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ABSTRACT 


The estimation of time varying spectra is a complicated one. The use of classical 
techniques coupled with the local stationarity assumption is met with only moderate 
success. Of the many time-frequency distribution functions used in the signal analysis, 
none present fully satisfactory spectra. The performance of the spectrogram, Instanta¬ 
neous Power Spectra (IPS) the Wigner-Ville distribution (WD) and various aspects of 
the Rihaczek distribution (RD) for a variety of signal nonstationarities are compared. 
WD has the most narrow main-lobes but suffers from spectral cross-terms. IPS, the real 
part of the RD consistently shows a bro jened main-lobe without cross-terms. The 
squared magnitude of the RD places sharp peaks along the crest of the main-lobe and 
is otherwise very similar to IPS. The imaginaiy part of the RD shows a sensitivity to 
discontinuous frequency changes i.e., frequen>. v shui.' eying. 
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I. INTRODUCTION 


The analysis of stationary spectra is a w ill-defined problem. Although from a the¬ 
oretical point of view a true signal spectrum can only be defined in terms of infinite du¬ 
ration data, spectral estimates resulting from the analysis of finite duration data have 
proven verj' useful. In the classical estimation problem, one assumes that 

• the data length is finite, and 

• the random process from which it originates is at least wide-sense stationary. 

Armed with these assumptions the behavior of spectra derived in this manner can be 
accurately predicted. The distortion incurred when analyzing a finite amount of data can 
be kept at a minimum given the specific estimation problem. But what happens if the 
signal is not at least wide sense stationarv? Signals commonly encountered in the real 
world are not stationaiy; they van.- in time. Either in amplitude or frequency content 
or possibly both, experimental data is rarely truly stationary. To better describe the 
variable random process, a time dependence must be included. All the theoretical results 
and even the practical application to finite duration data assume that ergodicity applies. 
Clearly a nonstationary signal is not an ergodic one, i.e., not one whose time average 
is equivalent to the mean realization. 

In the analysis of nonstationarj’ phenomena, there are a certain properties of the 
resulting spectrum which must be identical to the stationaiy' analog. These properties 
include an all-positive spectrum and zero energy in the spectrum when the signal is not 
present. Further constraints must be applied to the spectral behavior along the time 
dimension, a problem unique to ti e anaysis of nonstationary phenomena. There have 
been many attempts to adequately model the time-vaiy'ing behavior of nonstationary 
spectra. In general it appears that each technique has some advantages and disadvan¬ 
tages. Some appear better-suited to the analysis of a certain class of signals. There has 
yet to be found a completely satisfactoiy- description of the time dependent spectral es¬ 
timation problem. 
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II. SPECTRAL ANALYSIS OF SIGNALS WITH STATIONARY 

CHARACTERISTICS 

The Wiener-Khinchin theorem states that if j:(r) is a band limited, wide sense sta¬ 
tionary process, then the power spectral density (PSD), is related to the autocorrelation 
function (ACF), through a Fourier transform 

fxM = f” (1) 

‘'-oo 

where /?„(t) is the ACF. Using the Fourier inversion formula, 

= f ” (2) 

•'-OO 

and evaluating the correlation function at lag zero results in the average power in the 
process. 

where E denotes statistical averaging and 2T represents the duration of observation. 
Equations (1) and (2) describe the relationship between the time dependent data and the 
frequency dependent power spectrum. Equation (3) describes the relationship between 
a signal's temporal density and its spectral density. The integrand to the left in (3) re¬ 
presents the instantaneous power of the process. The integrand to the right represents 
the power as a function of frequency. Both can be considered power density functions 
and both are non-negative everywhere. Because the ACF is always conjugate symmetric, 
the power spectral density is always real and the average power in the process is always 
real. (Ref. 1,2) 
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Equation (1) can be rewritten, 


Pxxif) = -‘m £ 

7-^00 


2T 




”1 



2 

J 

-r 

- 


(4) 


which implies knowledge of the signal for all time. It is clear that in order to compute 
the true ACF or PSD, an infinite data set is required. Knowledge of an infinite number 
of realizations is also implied. These two theoretical constraints are impractical. Real¬ 
istically one must deduce the power spectrum from a single, finite duration realization. 
This chapter examines a variety of estimation techniques, noting the specific advantages 
and limitations inherent in each. 

A. CLASSICAL TECHNIQUES 

The computation of a power spectrum from one, finite set of data serves as an esti¬ 
mate of the true PSD. One conunon estimation technique is the periodogram. It can 
be computed from the data as 



(•T 

*'0 
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(5) 


an estimate similar in form to (4), except that limiting and statistical averaging oper¬ 
ations have been ignored. This estimate has the advantage of being both real and posi¬ 
tive. Examination of the mean and variance of the periodogram spectral estimate best 
describes its deviation from the true PSD, 


£[PfM] = 


sine I Tif - a)] PxJ^a)da, 


( 6 ) 


where it is apparent that the periodogram represents a smeared version of the true PSD. 
The smoothing along the frequency axis is caused by the finite observation interval. 
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Increasing the observation interval T narrows the main lobe of the sinc^ function, 
thereby minimizing the smoothing effect of the convolution in (6). The variance, on the 
other hand, is not so accommodating. [Ref 3] 

The variance of the periodogram spectral estimate for white Gaussian noise is a 
constant, with a standard deviation on the order of the mean. With so great a variance, 
the utility of the periodogram as described by (5) is questionable. One technique used 
to make this estimator more reliable is to compute an average periodogram. This has 
the effect of scaling down the variance by the number of terms in the average. In p»'ac- 
tice, since more than one realization is rarely available, this amounts to segmenting the 
data into shorter intervals, which causes a corresponding increase in the bias and loss 
of resolution of the estimate. 

One variation of the averaged periodogram scheme requires a data window. Look¬ 
ing at smaller segments of the data, the window is applied in an overlapped fashion and 
subsequent processing results in a series of periodograms that tend to be correlated and 
hence statistically dependent. Consequently, the actual reduction in variance will gen¬ 
erally be less than the number of terms in the average. Still another variation requires 
the data be prewhitened. This has the effect of reducing spectral bias, which is a prob¬ 
lem compounded by the averaging process. [Ref 3] 

Increasing the data length increases the resolution of the periodogram. An im¬ 
provement in variance can be realized if instead of a periodogram, one uses a 
Blackman-Tukey spectral estimator. This estimator is derived from a biased ACF esti¬ 
mate, 


j 4t)x\t -f r)dt, (7) 

where T is the duration of the data interval. Taking the Fourier transform of a win¬ 
dowed version of (7), 


A foO A „ , 

oo 

= J H'V- a)Pperia)da, 


( 8 ) 


leads to an improvement in variance. The type of window indicated by (8) should be 
capable of enhancing the spectral characteristics of interest [Ref 4]. An unavoidable 
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loss in fine spectral detail is the price paid to improve the variance of the spectral esti¬ 
mate. Here 2L and T are the durations of the lag window and the total data observation 
respectively. The variance of the Blackman-Tukey estimator is 

(9) 


y versus 

V'^ihMyplin 


( 10 ) 


for the periodogram. [Ref. 3] 

Classical techniques are reliable but have limited resolution and/or a poor variance. 
One assumes that the data is zero before and beyond the observation interval. This as¬ 
sumption causes the resulting spectral estimates to deviate from the true PSD. Using 
modern techniques, this type of error can be minimized. There are however, other limi¬ 
tations to consider. 

B. MODERN TECHNIQUES 

.Modern spectral estimation techniques rely on linear filter theoiy. Parametric 
modeling provides the foundation for many modern spectral estimation procedures. 
Rather than assuming the data to be zero beyond some arbitrary interval, parametric 
models assume statistical knowledge of the underlying process. The data is assumed to 
be composed of sinusoids in white noise. The model should accurately estimate filter 
coefficients for some linear filter whose output, when driven by white noise, is the 
available time data and hence has a spectrum exactly like that of the signal in question. 
The power spectral density for such an output is 

( 11 ) 

It is assumed that the transfer function, //(/), can be written as a rational polynomial 
and that the coefficients of the polynomial are the necessary filter coefficients with a\ the 
variance of the driving noise. Often the resonances of a particular random process are 
of interest. In this case an autoregressive (AR) model can be sufficient. 

AR modeling of spectra is the most popular modeling method. By solving a set of 
linear equations, accurate estimates of the parameters can be determined. Selection of 
an AR model results in an all-pole filter, hence the apparent high resolution. Selection 
of model order, p, determines the number of poles in the filter, which determines the 
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number of peaks in the spectral estimate. Once the coefficients are known, the 

AR spectral estimator becomes 

hM -T'*-:■ (12) 

1 + 

;=i 

The statistics of the AR spectral estimator are difficult to determine in closed form. 
Insight can be gained into its reliability by considering the following fundamental as¬ 
sumptions: 

• the process to be modeled is truly an AR process, 

p 

• the process contains — real sinusoids or p complex sinusoids, and 

• the coefficients {a(/')}f,i are accurate. 

The choice of model order obviously requires some knowledge of the signal's spectral 
content. The number of peaks in an AR spectrum, for complex-valued data, is equal to 
the model order, p. The spectrum, if the number of sinusoids truly present differs from 
p too much, can be unreliable. Furthermore, if the signal-to-noise ratio is poor, the es¬ 
timate can be of poor quality. [Ref 3] 

Other parametric models may be more appropriate to the process under investi¬ 
gation. Moving average (MA) models are used when the valleys or zeros of a spectrum 
are important. Unfortunately, MA models require the solution to a set of nonlinear 
equations. Although not impossible, it is much more dilficult to arrive at accurate pa¬ 
rameter estimates. Still another alternative is to use a combination of AR and MA 
modeling, resulting in what is called an ARMA model. 

The minimum variance spectral estimator (Capon's method) provides reasonable 
estimates without making any assumption about the data composition other than it 
originates from a wide-sense stationary random process. The classical periodogram 
spectral estimate can be interpreted as a bank of narrowband filters, all with identical 
passband characteristics. The minimum variance method is based on an adaptive 
passband characteristic with the advantage being the ability to adjust sidelobe levels and 
hence minimize spectral leakage. The resolution of this method lies somewhere between 
that of the periodogram and that enjoyed 'oy AR estimates. [Ref 3] 
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III. SPECTRAL ANALYSIS OF SIGNALS WITH DYNAMIC 

CHARACTERISTICS 


A. ADAPTING ESTABLISHED STATIONARY TECHNIQUES 
* In Chapter II, classical spectral estimators were found to be limited by the interval 

of observation, limiting the spectral resolution. A finite duration signal permits only a 
finite duration correlation estimate, with increasingly poor estimates away from the zero 
lag. Modern spectral methods have apparent higher resolution but are sensitive to the 
signai-to-noise ratio (SNR). All the estimators discussed thus far assume the data to 
be at least wide-sense stationary. The analysis of nonstationaiy phenomena complicates 
the esiimation problem. 

A wide-sense nonstationaiy process is one whose statistics or parameters var>’ with 
time [Ref. 5J. The actual nonstationarities of a process may include fluctuating power 
magnitude, changing frequency content or a combination thereof. There remains the 
problem of a suitable representation, one which appropriately displays the fluctuation 
in time and frequency of the instantaneous energy in the signal. To adapt stationary 
estimation techniques to the nonstationary case, the concept of local stationarity is iii- 
k troduced. A process is considered locally stationary if, over a given interval of time, the 

process appears to be stationary. Determination of the optimum interval depends upon 
the most rapidly fluctuating nonstationarity present in the process. If the process is 
observed overly long, the temporal fluctuations will be smeared. If the interval is un¬ 
necessarily short, spectral detail will be lost. A discussion of the more popular methods 
of nonstationary estimation and their limitations follows. 

1. Short-time Fourier Transform 

Short-time Fourier analysis is a method whereby the observed signal is seg¬ 
mented into a number of shorter intervals. Each segment is Fourier transformed and 
then magnitude squared. The resulting spectra are interpreted as cross sections of the 
true instantaneous spectrum. These cross sections are pieced together sequentially to 
form an estimate of the true time-varying power spectrum (Ref. 5, 6J. This type of 
spectrum is called a waterfall display by the signal processing community. Processing 
sequential sections of data results in a very crude estimate of the nonstationarities as a 
function of time. Furthermore, using longer segments increases spectral detail but tends 
to average or broaden the time-dependent fluctuations. An overly short interval will 
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tend to enhance the detection of the time transient behavior of the process at the ex¬ 
pense of spectral resolution. 

2. Spectrogram 

A modified version of the short-time Fourier transform estimation technique is 
the spectrogram It is related in that it uses a real, finite duration, sliding window 
(Ref. 7] centered at time t. The spectrogram 


hpearog.ir.‘)= ' 

*'-00 


(13) . 


computes a classical estimate of the spectrum for each point in time rather than for 
contiguous blocks of data (see for example. Figure 1). Similar to the periodogram, this 
spectral estimate is real and positive everywhere. The reliability of the spectrogram 
hinges on its ability to represent the signal's energy in the time-frequency plane. Refer¬ 
ring back to (3), any suitable time-frequency representation should reduce to 
when the frequency dependence is removed. Likewise, the representation should reduce 
to 1 X{J) I * when the time dependency is removed. Taking a look at the spectrogram's 
behavior. 




(a) amplitude (b) contour 

Figure 1. Behavior of the Spectrogram using an analytic sinusoid 
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mdf- 


|jr(T)| V(t -i)dx 


(14) 


and 


1*00 


(*00 


Spectrog 


mdi- 


\X\o)\^\lV{f-a)Vdo 


(15) 


we note that the estimate is smeared not only in frequency but in time axis as well, a 
result of the sliding window operation. Another indicator of the spectrogram's reliability 
can be found by removing both the time and the frequency dependency. The equality 
expressed in (3) for stationarv' phenomena should have a nonstationary counterpart. 
Therefore the volume under the spectrogram should equal the average energy in the 
signal: 


^oc 

^oo 

f*0O 


A 

^Spectrog.if^O^J^^ ” 



-00 ^ 

-oo^ 


^OO 


1 jf(T) 1 ^w^(t —t)dxdl. 


(16) 


The total signal energy requirement is satisfied only if the average energy of the window 
is equal to unity. The spectrogram represents a closer approximation of the instanta¬ 
neous energy changes relative to the short-time Fourier transform; however, temporal 
smearing is now unavoidable. [Ref 7, 8] 

3. Correlator/Matclied Filter 

The spectrogram estimator can be implemented as a bank of bandpass filters, 


Pcon^od) 



(17) 


X 1ZS« lb wviiw.iulCu V\iuil 


/y/0 = H< 


(18) 
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centered at each frequency of interest, and weighted according to some general lowpass 
characteristic, u'( —/) . The output of each bandpass filter is then magnitude squared, 
creating an instantaneous energy representation of the signal as a function of time and 
center frequency. Since the signal is known, the general passband characteristic can be 
replaced by a matched filter. Commonly used in a radar environment, each filter tests 
for round trip time delay and Doppler shift. The filters possess finite bandwidth and 
hence spectral resolution finer than dictated by the filters is impossible. Similarly, re¬ 
solution along the time axis is limited by the impulse response of the filters [Ref. 5). 

4. Autoregressive Modeling 

Autoregressive modeling can be adapted to fit slowly-varjing spectral charac¬ 
teristics. How slowly the frequency fluctuations must occur depends on the actual 
process in question. In general, as the signal rotates away from the pole of the AR filter, 
broadening of the peak results. As the correlation function is time-dependent, suitable 
accuracy in the AR coefficients requires an embedded time dependence. This particular 
modification is discussed in a later section. 

B. TIME-DEPENDENT SPECTRAL TECHNIQUES 

Adapting stationaiy techniques to the nonstationarv’ case is only marginally suc¬ 
cessful. There have been many attempts to describe the variation of signal energy, a 
function of both time and frequency, as a multivariable density. More appropriately, the 
simultaneous distribution of signal energy in time and frequency requires definition. 
Each technique has its merits, and each its own peculiarities. In 1966, a generalized 
phase-space distribution function was proposed [Ref 9], which can be used to derive 
many of the more popular time-frequency representations. This generalized distribution 
function is 


= J J J <I>(o, -h Y 


—oo OO oo 


• Y Vi) dt^ dr. (19) 


The advantage in using the generalized distribution (19) lies in the ability to define 
properties belonging to all representations derived in this manner. The distribution de¬ 
pends on the choice of O(o, t), referred to as the kernel of equation (19). An excellent 
presentation of the relationship between particular time-frequency distributions and (19) 
can be found in [Ref 10) and is summarized below. 
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1. The Running Spectrum 

a. Using Only the Past and Current Data 

In 1952 C.H. Page (Ref. 11] derived a time-frequency presentation arguing 
that one has knowledge of the signal up to and inci’ ding time t, but its future values are 
unknown. He defined a running transform, looKing backwards over all previous data 
as 


A7(/)=f' (20) 

OO 

where the superscript (-) indicates that the signal has been observed over the interval 
(- oo, t). By differentiating the squared magnitude of (20) with respect to time 

p-(/;r) = -|-|A7(/)l', (21) 

the running spectrum suggested by Page results. Substituting (20) into (21) and com¬ 
puting the partial derivative leads to the following alternate form 

= 2/?c[.v’(0A7(/)e^'"-^']. (22) 

M.H. Ackroyd has argued (Ref 5] that (22) is a finite duration approximation relative 
to the physical measurement of a true time-varying energy distribution. He suggests that 
the true time-var\ ing spectrum for a real-valued signal can be expressed as 

x{t)RelX<J)e’^’'^‘l = -^(0 f"" -^(v) cos(2nJ{t - T))dr, (23) 

OO 

the product of the response of a linear filter driven by the signal and tiie signal itself 
The implementation suggested by (23) is shown in Figure 2 (a). It requires an infinitely 
narrow filter with a noncausal impulse response. To deal with the causality issue, the 
impulse response is modified by a unit step function, 


j:(/) j Jr(T)«(r — t) cos{2nJ[i — t))</t = x(r) f x(t) cos{2nJ{t — x))dx 

oo oo 


= x(0Rc[A7(/)e''"-'''] 

=\p~m 


(24) 
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which is a scaled version of (22). Figure 2 (b) suggests a practical means of measuring 
a time-dependent spectrum. For ease of comparison. Figure 2 (c) shows the imple¬ 
mentation of a distribution discussed on pages 14-15. The distribution of a signal's en¬ 
ergy as proposed by Page is limited in spectral resolution; i.e., limited by the finite 
bandwidth of the filter. It represents an improvement over the spectrogram which was 
found to be smeared in both the time and frequency directions. 

Page's distribution can be generated from the generalized equation (19) us¬ 
ing the kernel function 

<I)(o,T) = e'''“'^', (25) 

which allows one to write (21) in yet another form. 



Equation (26) can be interpreted as the Fourier transform of an estimate of the true 
time-varying ACF, where 


R{t, t) = A- (/).v(r + t) for — CO < T < 0 
= x{i).x {t — t) for 0 < T < oo. 


(27) 


The behavior of this distribution can be seen in Figure 3 (a). As time increases, the 
amplitude of the signal continually increases and the true frequency location becomes 
more localized. 
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(a) amplitude (b) contour 


Figure 3. Behavior of Page''s distributions using an analytic sinusoid 


b. Using Only the Future and Current Data 

In 1967, M.J. Levin {Ref. 12] extended the concept of a running transform 

to include 


P CO 

Xt{f)=\ (28) 

where the superscript (+) indicates only those data values occurring at or later than time 
t are to be considered. The equivalent to (22) is 




(29) 


and the corresponding kernel function is 


(30) 


The behavior of this future term can be seen in Figure 4. Maximum frequency local¬ 
ization occurs early in the distribution, decreasing continually as time progresses, 
c. Using All of the Data 

Assigning equal weight to the past and future terms Levin defined the in¬ 
stantaneous power spectrum (IPS) as the average of the two running spectra. 


♦ 
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(a) amplitude 

Figure 4. Behavior of Levin's distribution for an analytic sinusoid 


= + A7(/)] 


Figure 2 (c) suggests a method whereby the IPS may be generated. IPS, like Page's 
distribution, requires a noncausal infinitely narrow bandwidth filter. Modifying the im¬ 
pulse response in order to create a realizable filter [Ref. 13, 14 [ causes both temporal 
and spectral smoothing, [Ref. 15, pp. 26-28]. The two terms of IPS can be interpreted 
as follows. The past term contains information of the energy and energy flow to create 
the signal up to time t. The future term contains the information about the energy and 
energy flow of the signal after time t [Ref. 6). 

2. Instantaneous Power Spectrum (IPS) 

IPS can be derived from the generalized time-frequency distribution using the 
kernel lunction: 


<l)(o, t) = cos(7toT). (32) 

This kernel, cos nor, can be formed by taking one half the sum of the kernels for the past 
and future running .spectra. Substituting (32) into (19) and simplifying gives 
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(a) amplitude 

Figure 5. Behavior of IPS for an analytic sinucoid 


(b) contour 


f*oo 


IPS(f,l) = j I 1 cos TtOT Jf(/, + — )x*(/j — ■— dt^ dx 

r-oo 

= yJ (jr(r) j:*(r - t) + x* {t) x{t + x))e~^^~^''dx 
= Re[x{t)X'{f)e~^^^^ 


(33) 


it is important to note that the terms m the integrand of (33) do not bear a direct re¬ 
lationship to the past and future spectra defined previously by (31). In this form, each 
term in the sum spans all the data and therefore contains contributions from both run¬ 
ning spectra (Ref. 15|. The ACF estimate as defined by IPS is 

RtPsiif = "o' (•’^(^) -"^ (^ ~ ^) + -*^*(0 + t) ) for — oo < T < oo. (34) 

By comparing Figure 3 and Figure 4, the behavior of IPS as shown in Figure 5, dem¬ 
onstrates an improvement in end-point resolution where the main ridge is most narrow 
at the center of the duration. 

3. Rihaczek Distribution (RO) 

Derived from physical consiuetaiions (Ref. 5,I6|, the complex energy or 
Rihaczek distribution (RD) is 
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( 35 ) 


Since the real part of a complex function is equal to the real part of the complex conju¬ 
gate of that function, it is obvious the IPS as defined in (31) is equivalent (33), the real 
part of (35). This relationship is depicted in i igure 2 (c). As yet, there is no satisfactory' 
interpretation of the imaginary part of the Rihaczek distribution (ImRD), although its 
computation leads to an increase in spectral localization over that of IPS for certain 
signals. A closer look at the behavior of IPS and RD can be found in Chapter IV. 

There exists a relationship between the spectrogram estimate and the RD. Re¬ 
writing (13) 


Spearog, 




/^oo 


jc(t)vv(t — t)e 


—OO 


f^oo 


.y(t)iv(t — 


—OO 

1^00 


.v(t)vv(t — l)e 




poc /»oo 


-OO -OO 

f^oo /•co 


X (ri)iv (tj — dx 

X\a)\V'{f- a)e>^^'^~'’'>‘da dx 
j:(t)A'' {a)e~^^^‘’'w{x — .'''IK (/* — dx 

RDJ^a, x) RDJf— a, x —t)da dx, 


(36) 


—OO —oc 


where RD{f,t) is defined by ^35). Equation (36) shows the spectrogram to be the 2-D 
convolution of the RD of the data with the RD of the window function. (Ref 6, 15) 

The complex energy distribution can be generated from (19) using the kernel 

function 


<I>(o, t) = 


(37) 


Substituting (37) into (36) results in 
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f^oo 


1^00 


RDm= 


x(f, + ^ )A:’(r, - -f ■‘"■•^•^i/i) dl^ dx 


4 / 4 / 

^-00*'-0C “0<> 

/^oo ^oo 


2 ' ' ' 2 
f^OO 


4t\ + Y - Y) 




dr 


(^oo 


x{i\ + Y ~ "f" i" dx 


( 38 ) 


—oo 
/*oo 


x{t)x (t — x)e 


= x[t)X{f)e 


■•(A.-P~'f‘ 


Equation (38) can also interpreted as the Fourier transform of an ACF estimate, where 

R{t, x) = x{t)x\i - x) for “ OO < T < OO. (39) 

Because RD is complex, this ACF estimate cannot be an even function of the shift var¬ 
iable, T. This suggests that it is the nonstationarities of a process which lead to an ACF 
which is partially odd. The behavior of this particular distribution is shown in Figues 5 
and 6. Figure 5 shows the behavior of IPS for an analytic sinusoid. IPS is equivalent 
to the real part of the RD. The imaginary part, shown in Figure 6 (a) and (b), demon¬ 
strates an improved sensitivity to rapid changes in signal energy as a function of fre¬ 
quency. The behavior of the RD is discussed in more detail in Chapter IV. 

4. Wigner-Ville Distribution (WD) 

Originally introduced by Wigner in a quantum mechanical context [Ref. 17] and 
extended to signal analysis applications by Ville [Ref. 18], the Wigner-Ville distribution 
(WD) is another valid representation of a signal's energy as a function of both time and 
frequency. The W’D can be generated from (19) using 

a)(u,T)= I (40) 

as the kernel function. Substituting (40) into (19) results in 
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(a) amplitude of imaginary part (b) contour of imaginary part 


Figure 6. Beiiavior of IniRD for an analytic sinusoid 


+ (41) 

‘'-oo 

The WD can be interpreted as the Fourier transform of an ACF estimate defined by 

t) = x'(t - Y ) x(r 4- Y ) for -oo<x<oo. (42) 

Since the WD is always real, this ACF estimate possesses even symmetry about the 
point of zero lag. An example of the WD is shown in Figure 7. 

5. Tiine-varj’ing Autoregressive Models 

An appropriate AR model for time-varying spectra contains an embedded time 
dependency. In so doing, the model would be able to track a spectral peak minimizing 
the effects of broadening. Rewriting equation (12) to reflect this time dependence, 

—*—-■ («) 

/=o 

where the coefficients are given by 
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Figure 7. . 


(a) amplitude 
Beliavior of WD fur an analytic sinusoid 


racoucNCi • ana 

(b)contour 



K 

k=0 

The {%}f=o represent the weights of a sum of ortho-normal basis functions. An appro¬ 
priate selection of ortho-normal basis set ideally uses some a priori knowledge of the 
spectrum under investigation. The time-varying AR model order is (/f -1- 1)P , requiring 
the estimation of {K -b 1) times more coefficients relative to the stationary analog. 
Whether one can accurately model the process with a manageable number of coefficients 
will depend on the particular process at hand. 
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IV. COMPARISON OF T-F DISTRIBUTIONS 


A. THEORETICAL RELATIONSHIPS 

Equation (36) shows that the spectrogram can be represented as the 2-D convo¬ 
lution of the RD of the signal with the RD of the window. A similar expression can be 
derived relating the spectrogram to the WD. Both RD and WD can be derived from the 
generalized distribution formula (19), as can many other t-f representations. It turns out 
that any particular distribution, Cif,t), whose kernel function satisfies 

<D(t),T)(I)(-u,T)= 1, (45) 

can be related to the spectrogram in the following way [Ref. 10], 


/•oo 


Spearog. 


m= 




C;t(a, t) Ch,(/’— a, T —i)do dr, 


oe —oo 


(46) 


where C, and Q are the generalized distribution functions of the signal and window re¬ 
spectively. The spectrogram itself, can be represented through the generalized equation 
using a kernel function 


<h(a 




w{[ + -y ) vt.' {l + 'Y 


'-v)t 


dt. 


(47) 


The right side of the equality in equation (47) is the ambiguity function (AF) of the 
window, w{i) [Ref 6]. For clarity, the AF is defined as [Ref 8] 


f OO ^ 

■^) = xil) (t -1- t)*; 

^-co 
r oo 

= - y) + Y 


(48) 
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So far, we have seen that a relationship exists between the spectrogram and certain 
other distribution functions (46), and that the spectrogram can be generated using a 
modified form of the AF as its kernel. In fact, all the distributions discussed thus far can 
be related through the AF. Rewriting (19) in a slightly difierent form, 

C{f,t) = [ f f 0(0, T)j:(r, + Y )x*(/, - y dt^ dr 

oo oo oo 

""1 J I 

r oo 

f~oo ^ 

where the generalized distribution is shown to be the 2-D convolution of the double 
Fourier transform of the kernel with the double Fourier transform of the complex con¬ 
jugate of the AF. Table 1 lists the distributions discussed thus far, along with their re¬ 
spective kernels. 


Table 1. DISTRIBUTIONS AND CORRESPONDING KERNEL FUNCTIONS 


Disiribulions 


Kernel - t) 

Specirogram 



Page 


iri 

Levin 

-fiArF 

IM 

1 

IPS 

Re [;c(0 X'{/) 

cos rruT 

RD 

x{t) X'(f) e-vi-A 


\VD 

IIIIIIIIIIII^IQQ^ 

1 


= f,,,[0(o,t)]’^;f;. 


-y2c(«+/t) 


dt) dr 


(49) 


B. GENERAL PROPERTIES 

Comparing Figure 1 and Figures 3 through 7, it is apparent that the particular 
kernel function has a tremendous influence on the particular properties of the resulting 
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distribution. From a spectral point of view, an appropriate t-f spectral representation, 
hence appropriate kernel Lnction, should ensure certain properties. These properties 
and corresponding kernel restrictions, if any, are listed in Table 2 below. 
[Ref 10, 15, 19, 20, 21] 

Of the many distributions discussed, none possess all the desired characteristics. The 
choice of a particular t-f representation depends on the application at hand. For prac¬ 
tical applications, the properties ascribed to the various distributions in Table 2 must 
be re-evaluated. When using windowed data, the resulting WD is referred to as the 
pseudo-Wigner distribution (PWD). Similarly modified versions of IPS and RD are in¬ 
dicated by the subscript (y), where y = jr(r) w(/) is the product of the data with some 
window function. Considering the linear operation of filtering, predictable distortion of 
the signal's true spectrum is encountered. In Table 3 a summaiy of the effects of four 
linear operations as they relate to time-dependent spectra is given. [Ref 10, 15, 19] 

C. RELATIVE PERFORMANCE 

In the previous section it was shown that the spectrogram is related to certain gen¬ 
eralized distributions through a 2-D convolution. It was further shown that any gener¬ 
alized distribution is the 2-D convolution of the double Fourier transform of the kernel 
function with the double Fourier transform of the complex conjugate of the ambiguity 
function (.AF). All these inter-relationships are interesting. By looking at the relative 
performance of certain t-f distributions, insight into the benefits and disadvantages 
characteristic of a specific representation are easily seen. The distributions compared in 
this section are: 

• Spectrogram 

• IPS, 

• ImRD,, the imaginary part of RD, 

• 2 linear, masnitude combinations of the real and imaginar\' parts of 

RD, 

• PWD, 

where the subscript (y) indicates the use of windowed data. 

1. Experimental Analysis 

Seven test signals, in the absence of noise, where considered. What follows is 
a brief description of the true time frequency behavior of the analytic signals and a list 
of the particular processing parameters of the digital implementations. Unless othenxise 
specified, all data is 128 points in duration, using a 1024 point FFT algorithm. Two 
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Table 2. GENERAL PROPERTIES OF T-F DISTRIBUTIONS 


Properties 

Kernel constraints 

Distributions 


xit) -* C(f,t) 


1 

2 


B 

5 

6 

Zero energ>‘ 

x{ta) = 0 -* C(f, to) = 0 

- --- .- . -. 


1 

B 

B 


B 

B 

X{fo) = 0 C{fo, t) = 0 


1 

B 

B 

fl 

B 

1 

Time shift 

“ ^o) t — (q) 

<I)(u. t) must be independent 
of absolute time and fre¬ 
quency 

1 

1 

S 

1 

\ 

1 

Freq. shift 

x{t) €^^"^ 0 ' -4 C(f-fo. t) 

Positive 

C(f,t) ^ 0 V/t 


X 

B 

B 

B 

B 

1 

Real 

C{f,t) = C{-f,t) 

0(o, t) = <I)’(-o,- t) 

B 

B 

B 

B 

B 

B 

Marginal in t 

f-C(f,t)df= |.v(0P 

-zo 

0) = 1 Vo 

1 

B 

B 

B 

B 


Marginal in f 

f-c{r,t)dt = iA-(/)p 

-oo 

(I)(0, t) = 1 Vt 

1 

B 

X 

B 

B 

fl 

Group delay 

/« t C{f,t)dt 

-00 Y' 

lAtOI* “ ‘ 

0(0,t) = 0 Vt 
and 

-^0(w,t)|,^o = 0 


1 


1 



Instantaneous^ 

frequency 

f-fC{f,t)df 

■u(OP 

0(0,0) = 0 Vo 
and 

0(i;, t) 1,^0= 0 


1 

1 

1 



Legend: 

1 = Spectrogram, 2 = Page, 3 = Levin, 

4 = /P5, S = RD, ,(>=WD 


1 signals only 
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Table 3. LINEAR OPERATIONS 



signal 

Fourier 

transform 

Distributions 

Ideal 

X(l) 


IPS,(f,t) 


WD,{f,t) 

lime shift 

x(l - to) 



RD.(f,i-k) 


modu¬ 

lation 


W-fo) 

iPS,{f-M 

RD,(f-M 


window¬ 

ing 

x(l) w(t) 

s 

* 

Re LRD,{f.t)*RD,m 

RD,{f,trRD,{f,i) 

WD,{f4riVDJU) 

filtering 

x(t) • h(t) 

m »'t/) 

Re lRD,{f,trRD,(f,t)-] 

RDM*RDM<t) 



difTerent Hamming window functions are used depending on the type of spectra ana¬ 
lyzed, stationar}', non-stationary, or a combination of both. Except for the spectrogram, 
all distributions have been smoothed in the time direction using a 5-cell box-car average, 
centered at the time of interest. The test signals used were: 

1. Single component., analytic sinusoid computed as where I < « < 128, using 

a 127 point Hamming window, 

2. Two component., analytic signal computed as where 

1 < « ^ 128, using a 127 point Hamming window, 

3. Singje component, analytic linearly chirped signal computed as 

, where 1 < « < 128, using a 55 point Hamming window, 

4. Two parallel., analytic., linearly chirped signals computed as 

where ! ^ m < 128, using a 55 point Hamming 

window, 

5. Sinele component, analytic, quadratically chirped signal computed as 

, where 1 ^ < 128, using a 55 point Hamming window, 

6. Multi-component signal comprised of a stationary’ sinusoid, a linearly chirped and 
a quadratically chirped component computed as 

, where 1<«<128, using a 55 

point Hamming window, 

7. Frequency shift keyed (FSK) signal computed as for 1<«<24 and 

56 <ii< 128, and for 24 ^ « < 56, using a 55 point Hamming window. 
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2. Highlights of the Analysis 

To judge the accuracy of component placement along the frequency axis, the 
performance of IPS, is compared to that of the spectrogram in Figure 8. Using a sta¬ 
tionary, single component, analytic signal both spectra have good end-point resolution. 



(a) IPS, 


(b) Spectrogram 


Figure 8. Amplitude plots of a stationary' analytic sinusoid 


Where the spectrogram presents a stationary spectral ridge with increasing amplitude 
near the center of the t-f plane, IPS, shows a fairly constant amplitude, stationary 
spectral ridge which is wider by comparison. Equation (36) describes the spectrogram 
as a 2-D convolution of the RD of the signal and the RD of the window. The apparent 
superior resolution of the spectrogram over IPS,, the real part of the RD, indicates that 
the imaginary part of the RD contains some important spectral information. The im¬ 
aginary part is formed as 


ImRD(f,i) 


X 

2 


^OO 


(.v(/) X* {i — t) — x{t + 


(50) 


Using the kernel function 


tI>(o, t) sin TTOT, (51) 

equation (50) can be derived from (19). It is obvious that the kernel for the RD is the 
sum of the kernels for IPS and ImRD. As to the behavior of IntRD,, a windowed version 
of (50), demonstrates an improved sensitivity to spectral change relative to PWD and 
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(a) Amplitude (b) Contour 


Figure 9. ImRDy for a statiunar}', analytic sinusoid 

IPS^ whicli is shown in Figure 9. In particular, using a single-component, linearly 
chirped signal, the zero crossings of ImRDy occur near the true frequency location. In 
Figure 10, the relative performance between PWD, IPSy and ImRDy for a single¬ 
component, linearly chirped, analytic signal can be compared. One can see that ImRDy 
has an improvement in resolution over that of IPS, comparable to that achieved by 
PWD. Where both IPSy and PWD rely on the sharpness of the spectral ridge for resol¬ 
ution, ImRDy relies on the zero crossings. 

The ImRDy has large values near the time of frequency change, with the ampli¬ 
tude approaching zero otherwise. This behavior is clearly demonstrated in Figure 9 
where a single component stationary signal is shown to have a nonzero imaginary 
spectra, greatest in amplitude near the beginning and end of the life of the signal. In 
Figure 8 (a), IPSy presents a spectral ridge which, although rounded at the endpoints, 
grows increasingly more narrow toward the center of time. By forming a linear combi¬ 
nation of IPSy and ImRDy, an overall improvement in spectral resolution for some 
signals can be achieved. One such linear combination is 

I RDy{f,t )!' = I IPS^{f,i) 1 2 + I ImRDy{f,t) I ^ (52) 

Another combination can be formed taking the difference of magnitudes. 


\RDy[f,t)t=\IPS^if,t)\'^ - \ImRDy{f,t)\\ 


( 53 ) 











ncojocr • mis 

(c) Contour plot of ImRD, 


Figure 10. Contour plots of a single-component linear chirp 

Equation (52) results in a nonnegative spectrum. Equation (53), unfortunately, can have 
negative values. 

For stationary data, represents an improvement in the end-point resol¬ 

ution relative to IPSy as can be seen in Figure 11. When two stationary components 
are present IPS,, ImPD, and both linear combinations produce modulated spectra. The 
modulation effect is related to the difference frequency. PWD produces a spectrum 
having cross-terms midway, between the true components. The cross-terms are also re¬ 
lated to the difference frequency. The relative behavior of IPSy, and PWD for 
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Figure 11. Contour plots of a single-coinponent, stationary analytic sinusoid 

a two-component analytic signal is shown in Figure 12. Using 128 point data sets at a 
sampling frequency of 128, classical analysis predicts component resolution where the 
separation is at least 0.89 Uz (Ref 4). IPS^ is capable of resolving two components 
separated by at least 1.25 Hz. |i?/),|i and PWD are able to resolve two components 
separated by at least 0,6 llz, an improvement over classical analysis. Figure 13 is a 
graph depicting the relative location of the spectral peaks as a function of frequency and 
time. IPSy, after an initial settling, consistently places the spectral peaks demonstrating 
a bias which is not symmetric. similar to IPS^, consistently places the spectral 

peaks but with a symmetric bias. PWD succeeds in making the nearest approximation 
to the true component locations demonstrating a bias which is not symmetric. The 
placement of these spectral peaks does not appear to settle at one location as appears 
to be the case for the Rihaezek-derived distributions; however in a -^ean-squared error 
comparison, PWD appears to be superior. 

The ImRD, is also capable of resolving two closely-spaced, narrow-band, sta¬ 
tionary components. Instead of searching for spectral peaks, ImRO, characteristically 
detects the zero crossings which yield information in this spectrum. In Figure 14, the 
behavior of ImRDy using two closely-spaced stationary sinusoids can be seen. 

To study the behavior of any spectral estimator of nonstationary phenomena 
we begin by considering a single-component, linearly chirped, analytic test signal. The 
relative perlormance ol IPS^, ImRDy and PWD was discussed previously, see also 
Figure 10. The behavior of the two linear magnitude combinations can be seen in f'ig- 
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(c) Amplitude plot of | RDJi 

Figure 12. Amplitude plots of a two compouent stationary, analytic sinusoid 

urc 15. Forming the dificrence, creates a spectral ridge comparable in width 

to IPSy but with a better defined peak. Detection using \RDy\i should be thus be im¬ 
proved. Forming the sum, | RD^V completely resolves IPS into its past and future terms, 
the distributions defined by Page and Levin. Component location using 1 requires 
detection of the minimum occurring between the two resolved ridges, making detection 
using \RD,\^ questionable. In the presence of noise, the separation IPS, into compo¬ 
nent parts will lead to difficulty in interpretation. The ability of IPSy to properly locate 
the instantaneous frequency for a linear chirped signal is shown in Figure 16. The lo- 
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Figure 13. Gmph depicting accuracy of stationary component placement 

cation of the peaks in /AS, coincide with the true instantaneous spectral locations for a 
slow chirp. Eatly in the spectrum, when the future term in dominant, IPS, tends to place 
the instantaneous frequency higher than truth. On the other hand, late in tnc spectra 
when the past term is dominant, IPS, lends to place the instantaneous frequency lower 
than truth. Doubling the chirp rate results in a greater frequency ambiguity as seen by 
tlie two terms which compose IPS,. Looking at rigurc 16, this is demonstrated by the 
apparent random placement by IPS, of the spectral peaks, neither the future nor past 
term seems to be favored as the ma,\inu’m peak location. 

Next, a test signal composed of two, parallel, linearly chirped, analytic signals 
is considered. 1 lie spectra resulting from IPS, and I’WD is shown in Figure 17. Similar 



Figure 14, IiiiRD, for a tuo-compoiient stationary, analytic sinu.. .d 
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(a) Amplitude plot of 1 y?Z)^ I i (b) Amplitude plot of | RDy \ ^ 

Figure 15. Relative perfunnance fur a single-coinpoiieiit linear chirp 

to the stationary, two-componcnt case, modulation and cross-terms fluctuating at the 
difference frequency are present in IPSy and PWD, respectively. Neither of the two lin¬ 
ear combinations of magnitudes of the component parts of the Rihaezek. distribution, 
was able to improve on the resolution achieved by IPSy. Moving from the peaks to the 
valleys, ImRDy can be seen to discriminate between the two chirps as a function of the 
zero crossings. This is demonstrated in Figure 18 where a pattern of zero crossings al¬ 
lows the eye to discern the two components. 



Figure 16. Graph depicting accuracy of IPS in locating the instantaneous frequency 
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(a) Amplitude plot of IPS, (b) Amplitude plot of PWD 

Figure 17. Beliavior of IPS, and PWD for 2 parallel linear cliirps 


Hxamination of the spectra for more complicated signal suggests a relative 
ranking in terms of resolution among the dilfcrent power distributions discussed thus far. 
1 he test signal is composed of a high frequency stationary component and two chirped 
components, one a linear chirp and the other a quadratic chirp. Considering first those 
spectra which estimate frequency location as the point of ma.ximum power, PWD 
produces the most narrow ridges; sec Figure 19 (a). It sullers from poor end-point re¬ 
solution and spectral cross-terms. seen in Figure 19 (b), produces well-defined 

periodic peaks along the instantaneous frequency path of the nonstationary compo¬ 
nents. The stationary component ridge is broadened, with the modulation effect most 


Figure 18. 



ImRDy of tuo parallel linear chirp.s 
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(c) ImRD, 

Figure 19. A coinbinatiun of stationary and nonstationary spectral components 

apparent near the end-points. !/?£>, P would be dilHcuIt to interpret as it resolves the 
past and future terms for the nonstationary components. Looking for the zeros, ImlU), 
accurately describes the location of the nonstationary components, but provides little 
information during parts of the stationary signals e.xistence. Again, Figure 19 (c) re¬ 
quires pattern recognition to be able to discern the individual, dynamic components. 

Considering a pulsed spectra, such as in FSK modulation, \RDy\^ presents a 
narrow ridge with good resolution throughout the duration of each pulse. 1 he ridge 
width for each pulse in the I’WD spectral description is dependent upon the pulse du- 
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(a) (b)PWD 

Tigure 20. Contour plots for a complex anlaytic FSK signal 

ration, 'l lic longest pulse showing the most narrow ridge, one which is slightly inore 
narrow than what is found with \RDJ^. These two spectra are shown in Figure 20. 
rWD has a slow build up and decay at the ends of each pulse. There appears to be a 
ttadc-oir between the width ol'a spectral lidgc aitd cnd-p('int resolution. 

3. Test Case Results 

Figures of the spectral dist. ibiitions resulting from the six t-f representations are 
contained in this section, 'fhey arc oidcied by the type of test signal under analysis. 
a. Single-component, Analytic Sinusoid 

Using the spectiogram as the reference, the resolution ability of five addi¬ 
tional t-f distributions is compared in Figure 21 - Figure 24. I’WD presents a well- 
defined spectral ridge near mid-plane, but sulfers from a sluggish built-up and decay. //’S, 
presents a wider main-lobe relative to I’WD which is compensated somewhat by an in¬ 
crease in end-point resolution. provides the best end-point resolution, main¬ 

taining constant amplitude throughout the plane, |/?/.), I i proves to be more sluggish 
than PWD and IniRD, demonstrates an improved response near the time of spectral 
change. 
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Figure 22. Test sign.il I: contour plots for Spectrogram, IPSy and PWD 
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(c) ImRDy 

Figure 23. Test signal 1: ainpliitide plots for 17?Z>ylS and ImRDy 
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b. Two-component, Analytic Signal 

The spectrogram and all five t-f distributions shown in Figure 25 through 
Figure 28 display a distorted spectrum. This distortion exists for a minimal time in the 
spectrogram estimate. In the case of IPS, and PWD, the distortion is related to the 
dilTerence frequency. IPSy shows modulation of each component; PWD contains addi¬ 
tional peaks oscillating at the difference frequency. The modulation effect characteristic 
of Rihaczek-derived distributions makes the spectra using 1 RD, \ ’ and ImRDy extremely 
difficult to interpret. suffers also from the modulation effect; however, this 

magnitude combination enhances resolution relative to IPSy, making it easier to detect 
the presence of two components. 

c. Single-component, Analytic, Linearly Chirped Signal 

The inadequacy of assuming local stationarity is clearly demonstrated in 
Figure 29 (a) and Figure 30 (a) for the spectrogram, where the slope of the instantane¬ 
ous frequency line is distorted and broadened near the end-points. Both IPSy and PWD 
shown in Figure 29 (b) and (c), and Figure 30 (b) and (c) maintain a better approxi¬ 
mation to the instantaneous frequency slope neat the end-points. Characteristically, 
IPSy presents a broadened spectral ridge whereas PWD deca> s to zei o at the start and 
stop of the chirp. Looking at the spectrum created by I/?£>,.p. Figure 31 (a) and 
Figure 32 (a), the future and past terms which make up the RD are clearly visible.. Al¬ 
though 1 RDy I ^ is an all-positive distribution possessing many desirable properties, this 
characteristic resolution quality makes it difficult to interpret more complicated spectra. 
The two remaining Rihaczek-related distributions represent improvements over the 
ability of IPSy to pinpoint the instantaneous frequency as a function of time. | RDy \ i 
in Figure 31 (b) and in Figure 32 (b) shows a spectral ridge more narrow than that for 
PWD. ImRDy. while accurately locating the instantaneous frequency, requires the de¬ 
tection of zero crossings. For particular detection schemes information thus provided 
may be appropriate. 

d. Two Parallel, Analytic, Linearly Chirped Signals 

For this test signal, the spectrogram spectral estimates shown in Figure 33 
(a) and Figure 34 (a) is unacceptable. The comments made previously concerning the 
distortion in the spectra when when two closely spaced, parallel, stationary components 
can also be applied to this nonstationar)’ case (see Figure 33 - Figure 36). In the case 
of PWD, not only are cross-terms oscillating between the true components present but 
the spectral ridges themselves show the effects of modulation. |/?Z)yP and ImRDy do 
not show promise as estimation tools for close!}-spaced frequency components. It is 
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interesting to note however, that resolve the past and future terms for each 

component similar to the behavior shown in Figure 31 (a), \RDy\ \ shows very sharp 
peaks along the modulated instantaneous frequency lines of the two linear chirps. Be¬ 
cause these peaks are the largest peaks in the plane, | RDy\ \ from a practical view point, 
appears to be more suited toward the anaKsis of this type ol signal liian PWD whcic uic 
cross terms periodically represent the dominant peak. 
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Figure 26. Test signal 2; contour plots for Spectrogram, IPSy and PWD 
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(c) ImRD, 

Figure 27. Test signal 2: ainplitutle plots for | RDy\^, \ RDy\i and linRDy 
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(c) ImRDy 
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Figure 28. Test signal 2: contour plots for i/?Z)ylS and ImRDy 
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(c) I’WD 

29, Test signal 3: ainplitude plots for Spectrogram 

























(c) ImRL\ 

Figure 31. Test signal 3: amplitude plots for \RDy\i and ImRDy 
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Figure 32. Test signal 3: contour plots for liRD^lS \RDj,\i and ItiiRDy 
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Figure 34. Test sigii.il 4: contour plots for Spectrogram, IPSy and PWD 
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(b) \RD,\i 


(c) ImRD, 

riguie 36. Test signal 4: contour plots for 1 RD^W \ RDy\i and linRDy 


'■-2 







c. Single-component, Analytic, Quadratically Chirped Signal 

The ability of the spectrogram and these 5 t-f distributions to accurately 
display more rapid spectral dynamics can be compared in Figure 37 to Figure 40. Not 
surprisingly, the spectrogram presents a poor estimate. IPSy demonstrates a spectral 
ridge which tends to broaden as a function of time, making it difficult to ascertain the 
actual instantaneous frequency curve. PWD, although zero at the end-points, tracks the 
chirp closely presenting a ridge along the line of instantaneous frequency as narrow as 
that found for the linear chirp. appears to be provide the most narrow ridge. 

Both PWD and l/?Z)^|i show amplitude modulation along the peak of the curve. As to 
the remaining spectra, 1 RD^^ and ImRDy behave in a manner similar to the case of the 
linear chirp and thus do n^ appear particularly suited to this class of signal. 

/. Multi-component Analytic Signal 

How the various spectral estimation techniques perform when confronted 
with a mixture of stationarx' and nonstationaiy dynamics is demonstrated in Figure 41 
through Figure 44. This test signal is composed a high frequency stationary component 
and two chirped components, a linear chirp and a quadratic chirp. Their relative per¬ 
formance suggests a ranking in terms of desirability as an estimator of continuously 
changing spectral components. The spectrogram is predictably, the worst, followed by 
and ImRDy in ascending order of desirability. Only IPSy and PWD 

present spectra which closely resemble the true signal components. All three show a 
broadened ridge for the stationaiy component relative to the single stationaiy compo¬ 
nent case examined previously. This results from using a shorter window. As expected, 
I RDy I i sharpens the modulation peaks found in its own spectrum and that of JPSy. A 
comparison of | RDy | i and PWD is also quite characteristic. The price paid for in¬ 
creased end-point resolution and elimination of the spectral cross-terms is a broadening 
of the spectral ridge and the appearance of modulation along the crest of the spectral 
ridges. 

g. Complex FSK signal 

In contrast to the continuously-varying frequency dynamics of the previous 
test case, a veiy dilTerent order of desirability is suggested in the case of pulsed spectral 
dynamics. Comparison of Figure 45 - Figure 48 shows \RDy\^ to possess superior 
end-point resolution ability coupled with a rapid build-up and decay of the spectral ridge 
reiatise to ihe irue pulse dsnamics. IPSy and produce siiiiilai spectra; however, 

the end-point build up is slower. ImRDy, with its heightened sensitivity to detect spectral 
change appears to be an excellent indicator of the time and location of frequency 


5.3 







0 fWflUDCI MIS ■ 


(c) PWD 

Figure 37. Test signal 5: amplitude plots for Spectrogram, IPSy and PWD 

change. PWD presents a sluggish build-up and decay, with the wridth of each spectral 
ridge dependent on the duration of the pulse. Furthermore, cross-terms can be seen 
between the two higher frequency pulses. 
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(c) PWD 

Figure 38. Test signal 5: contour plots for Spectrogram, IPSy and PV\ D 
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Figure 40. Test signal 5: contour plots for | RDy \ \ | RD^ 11 and ImRDy 
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(c) PWD 

Figure 41. Test signal 6: amplitude plots for Spectrogram, IPSy and PVVD 
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(a) Spectrogram 
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(c) ImRD, 

Figure 43. Test signal 6; ainplitiule plots for \RDy\^, |y?Z)^lL and ImRDy 
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Figure 44. Test signal 6: contour plots for | RD^ | RDy\i and IniRDy 
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(c) PWD 

Figure 45. Test signal 7: atnplitude plots for Spectrogram, IPS^ and PWD 
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(c) ImRDy 

Figure 48. Test signal 7: contour plots for | 1 RDy\i and IinRD^ 
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D. ALTERNATE METHODS OF COMPUTING IPS 

Equation (33). the defining equation for IPS, can be rewritten as the Fourier trans¬ 
form of an ACF estimate (34). namely 


f*oo 


IPS{f,t) 


± 

2 


(x{t)x (r — t) + X (r)j£’(r + dx. 


/'oo 


(54) 


IPS for finite-duration discrete signals is 


IPS{d,n)=^^‘ Yj {x{n)x\n-k) + x\n)x{n + k))e~^^'', (55) 

/C=—OO 

where the signal sequence j:(rt) is finite and zero outside the known samples, and AT is 
a constant. Equation (55) can be expressed as (Ref. 15] 

lPS{d,n) - I x{n) i' = I X{d) I' - I X,{e) I(56) 


where 


OO 

we) = Y 

r=~oo 

r=—OO 


(57) 


For use later in the derivation, let the Fourier transform of the point of interest be 

m=A«)r'‘" .... 

iwjp-uwr. 

By neglecting U(.v)p in (56) a modified version of IPS is defined. The behavior of this 
modified IPS for a single-component, anahtic sinusoid is shown in Figure 49 (a). 
Comparing this with an unwindowed version of IPS for the identical data sequence in 
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(a) Modified IPS (b) IPS 

Figure 49. IPS 

Figure 49 (b), shows the two methods of computing IPS to be very similar. The greatest 
difference can be seen when comparing maximum amplitudes. Modified IPS is on the 
order of 10"’ and IPS is on the order of 10. Both methods were computed using a five- 
cell-box-car averaging procedure along the t axis. 

The question arises concerning the implementation of modified IPS using a window 
function. One method of windowing can be written as 


oo 

IPS{d,n) = ^ — A) -t- X*(n)x{n + k))w{k) (59) 


OO 


Applying the definitions in (57), 


= -^{D{d){x\e)yvie)) + z)*(O)(^(0)*if(0))} 


AT 

2 


OO ^ 

I D{0) ^ X^ipWiO -p) -t- D\0) £ XipWiO -p)l, 


CO 


(60) 


where * denotes convolution in the frequency variable. Moving everything inside the 
$ 

summation sign and substituting Xl(p) + D'ip) for X'(p) and A'*(0) — X'XO) for D'{0) gives 
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oo 

vn 


IPS{e,n) = -^ ) (A';(p) + D(p))D{d)lVid-p) 


LU 

p=«-oo 


y (Ai>) - XmUwno-p) 


P=—OO 


) [o(())x;w + mo'(f)]ine-p) 


OO 


( 61 ) 


P=—30 


+ 


AT 


oo 

^ [A'*(0)A'(p) - ^;(0)A'(;;)]in0-/.). 


p=-oo 


A straight-forward relationship between IPS and the spectral contribution of any single 
point in the data set is given by equation (56). Equation (61) defines an analogous re¬ 
lationship between (57) and (59). Unfortunately, this relationship is not so straight¬ 
forward and requires additional analysis to determine the benefits, if any, to be gained 
from processing data in this manner. 
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V. RECOMMENDATIONS AND CONCLUSIONS 


The comparative behavior of the spectrogram, IPS, WD and three novel t-f distrib¬ 
utions was explored in Chapter IV. Of the three, ImHD^ proves to be particularly sen¬ 
sitive to discontinuous changes in frequency. This characteristic behavior may be useful 
in certain detection applications. When confronted with more complicated spectra, one 
containing closely-spaced stationary^ components or continuously-varying nonstationary 
components, IniRD^ does not appear particularly useful. Similarly provides an 

improvement in end-point resolution when used for the estimation of stationary or 
pulsed spectra. Taking the square root of reiuits in a t-f distribution, denoted 
by II, which satisfies many of the desirable properties listed in Table 2. In 
particular |/?£)(/» j : 

1. Satisfies both zero energy requirements 

2. Obeys the time and frequency shift properties 

3. Is positive and real for all time and frequencies. 

How well i RD, | - can track a rapidly fluctuating pulsed signal, similar to something 
found in frequency hopped communications, is an area worth investigating. 

For detection of continuously changing spectral dynatnics. \RD^\i appears to be a 
viable processing scheme, however the performance of this estimator in noise needs still 
to be examined. Using \RDJl as an estimator is not v,'ithout problems. Unlike 
|/?/)|, taking the square root of |/J/)li does not produce an estimator which satisfies 
the marginal requiicnients. It does however comply with the zero energy and 
time frequenq- shift properties desirable for time-dependent estimation problems. Used 
with stationary spectra \RDy\i appears to improve the resolution of closely spaced 
components beyond that currently achieved by classical methods. For stationary 
spectra, 1 RDJi is a biased estimator, i.s., the true spectral peaks occur at a given fixed 
offset from their true frequency location. PWD also provide.s a resolution improvement; 
however the estimate it provides never settles down to one frequency location. For short 
duration data may prove superior in the detection of multiple stationary com¬ 

ponents. 

In addition to the experimental results presented in Chapter IV, /F5y demonstrates 
a 3-dB noise advantage relative to PWD (Ref. 15J. Coupled with superior end-point 
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resolution IPS may be the desirable method in practical analysis problems. All the re¬ 
sults presented were derived from noise-free data sequences using digital implementa¬ 
tions (see Appendix A). Because ImRDy, \ and | /?Z),|i all show advantages when 
applied to the appropriate signal type, analysis of their noise performance is an open 
issue. As they are derived from the RD, as is IPS, it is likely that they enjoy a similar 
robustness. 

Looking to practical applications of IPS, a brief discussion of performance in a 
multi-sensor environment can be found in Appendix B. The initial results look promis¬ 
ing, but more extensive research needs to be conducted. A second practical application 
scheme involves the use of a cumulant or third-order moment. Typically associated with 
the sonar environment, this scheme seeks to take advantage of the fact that the odd 
moments of a zero-mean, gaussian noise process are identically zero. An initial investi¬ 
gation using IPS to compute cumulants can be found in Appendix C. 

In general, WD produces ver\' narrow spectral ridges but suffers form poor end¬ 
point resolution and spectral artifacts. IPS provides an improvement over these short¬ 
comings at the cost of spectral broadening. Recently appearing in the literature is 
another t-f distribution [Ref. 10, 22j. Defined by H.l. Choi and W.J. Williams, this t-f 
distribution minimizes the effects of the spectral cross-terms. Closer examination of the 
resultant spectra, which uses the kernel function 

2 2 

<l)(y. x) — e 2t7 , 

where is a constant, shows that the spectral ridges are broadened, similar in behavior 
to IPS. In classical estimation, i.e., the periodogram, one maximizes spectral resolution 
using nothing but the raw finite data set. The price paid is a large, slow roll-off sidelobe 
structure that can mask a true component. In time-frequency distributions, i.e., in the 
generalized phase-space equation (19), one maximizes spectral detail using a constant 
kernel of unit amplitude. The price paid is poor end-point resolution and spectral 
artifacts across time and frequency. In the classical analog, using any other window 
results in improved sidelobe behavior at a corresponding loss in detail along the fre¬ 
quency axis. IPSy, the Rihaczek derived distributions and the distribution suggested by 
Choi and Willirnas all improve or eliminate the disturbing spectral cross terms charac¬ 
teristic of WD. The price paid is loss in spectral detail. 
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APPENDIX A. COMPUTER CODE 


To conserve space only the read file and four of the si.. Fortran codes have been 
included. Data uas generated and read by the basic programs as required. The data 
generation file is not included here. The read file allows easy change of processing pa¬ 
rameters. Further, IPS, and ImRD, are generated by code that difiers in a minus sign 
in line code. IPS, requires a plus sign; IniRD, requires a minus sign. Similarly i 
and are simply related. The location of the sign change is indicated in the re¬ 

spective algorithms. Graphics are produced using DISSPLA. 

1. Parameter File 


1 1 031 039 08 

HAMMING WINDOW ($' 

B 

IPS:$' 

FS=1/128, 40 POINTS OF DATA$' 
5 CELL TIME SM00THING$' 


MODE PLTR BWLEN EWLEN WINC 

WTYPE 

CONTR 

TTL 

SIGNAL 

TSMTH 
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2. Spectrogram 


c** this fortran file computes the spectrogram of *** 

C** A DATA SEQUENCE *** 

* 
* 

INPUT DATA SEQUENCE IS READ USING FILEDEF 4» AS THE COMPLEX * 
ARRAY DATA(L). * 

* 

L IS THE LENGTH OF THE DATA SEQUENCE AND IS ADJUSTED FROM THE * 
PARAMETER STATEMENT. L MUST NOT EXCEED 128. * 

if 

ANALYSIS PARAMETERS ARE READ USING FILEDEF 41. THE PARAMETERS * 
ARE: * 

ARGUMENT TYPE ALLOWED VAULUES * 

* 

MODE - II 1 PLOT 0 TO PI * 

2 PLOT PI TO PI * 

* 

PLTR - II 0 SHERPA LASER PRINTER 


1 IMB79 GPJ^PHICS TERMINAL 


BWLEN- 13 3 DIGIT INITIAL WINDOW LENGTH, * 

MUST BE AN ODD INTEGER * 

EWLEN- 13 3 DIGIT FINAL WINDOW LENGTH * 

WINC - 12 2 DIGIT WINDOW INCREMENT, MUST * 

BE AN EVEN INTEGER * 

* 

WTYPE- A19 19 CHARACTER STRING USED IN THE * 

PLOT HEADER DISCRIBING THE * 

WINDOW USED. THE CURRENT * 

WINDOW LENGTH IS AUTOMATICALLY * 
INCLUDED * 

Vf 

CONTR- A1 1 CHARACTER STRING INDICATING * 

TYPE OF PLOT DESIRED * 

* 

A AMPLITUDE PLOT ONLY * 

C CONTOUR PLOT ONLY * 

B BOTH AMPLITUDE AND CONTOUR * 

* 

TTL - A43 43 CHARACTER STRING USED IN THE * 

HEADING WHICH DESCRIBES THE * 

ALGORITHM AND THE CLASS OF * 

SIGNAL USED * 

* 

SIGNAL- A43 43 CHARACTER STRING DESCRIBING * 

TEST SIGNAL * 

* 

OUT REAL * 

OUTPUT ARRAY OF DIMENSION 512 BY L * 

V* 

COEF COMPLEX 

ARRAY OF DATA AFTER IT HAS BEEN WEIGHTED WITH A * 

SLIDING WINDOW FUNCTION. * 
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* 

FT COMPLEX * 

ARRAY OF THE 1024 POINT TRANSFORM COEF * 

* 

Z INTEGER * 

LENGTH OF CURRENT WINDOW * 

* 

M INTEGER * 

MID-POINT OF THE CURRENT WINDOW * 

* 

AMAX REAL * 

MAXIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS * 

* 

AMIN REAL * 

MINIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS * 

* 


* 


PARAMETER(L= 32) 

INTEGER I,J,N,M,MODE,Z,BWLEN,EWLEN,WINC,PLTR 
REAL OUT(512,L),AMAX,AMIN 
CHARACTER WTYPE*19,TTL*43,SIGNAL*43,CONTR*l 
COMPLEX DATA(L),FT(1024),COEF(1024) 

CALL EXCMSCFILEDEF 4 DISK TEST IN A (PERM') 
CALL EXCMSCFILEDEF 41 DISK PARAM IN A (PERM') 


--- READ IN PARAMETER LIST . 

READ(41,400)MODE,PLTR,BWLEN,EWLEN,WINC,WTYPE,CONTR,TTL, 

+ SIGNAL 

400 FORMAT (1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X,A43/ 
+ 1X,A43) 


TEST TO ENSURE WINDOW LENGTH IS APPROPRIATE . 

I=L-1 

IF ((BWLEN .GT. I) .OR. (EWLEN .GT. I)) THEN 
WRITE(*,69) 

GO TO 99 
ENDIF 

I=M0D(BWLEN,2) 

K=MOD(WINC,2) 

IF (I .EQ. 0) THEN 

IF (K .EQ. 1) THEN 
WRITE(*,68) 

GO TO 99 

ELSE 

WRITE(*,67) 

GO TO 99 
ENDIF 

ENDIF 

69 FORMAT (IX,'WINDOW LENGTH EXCEEDS LENGTH OF THE DATA') 

68 FORMAT (IX,'WINDOW INCREMENT MUST BE EVEN') 

67 FORMAT (IX,'INITIAL WINDOW LENGTH MUST BE ODD') 

C. 
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c 

C.PLOTTING DEVICE CALL 

IF (PT-TR .EQ. 0) THEN 
CALL COMPRS 
ELSE 

CALL IEM79 
ENDIF 


PI=4*ATANC1. ) 

READ( 4, *) ( Dii TA( I), 1=1, L) 


DO 111 Z=BWLEN,EWLEN,WINC 

M=(Z-l)/2 

CALL ANGLE(0.0) 

AMAX=0 
AMIN=AMAX 
DC 10 J=1,L 
DO 20 N=-M,M 

IF ( (d+N) .6E. 1) .AND. ((I+N) . LE. L) ) THEN 
COEF(M+N+l)=DATA(I+N) 

+ -'KO. 54+0. 46*C0S(2*PI*N/(2-^M))) 

ELSE 

C0EF(M+N+1)=(0. ,0. ) 

ENDIF 

20 CONTINUE 

DO 30 N--2*M+2,J.024 
COEF(N)=(0. ,0. ) 

30 CONTINUE 

CALL FFT(1024,COEF,FT) 

IF ( MODE.EQ. 2) THEN 
DO 40 N=l,513,2 

OUf(INT((N+1)/2+255),I)=ABS(FT(N))**2/(2*M+1) 
IF (0UT(INT((N+n/2+255),I) • GT. AMAX) THEN 
AMAX=OUT(INT((N+1'/2+255),I) 

ENL^F 

1? (0 ■r(INT((N+l)/2+255),I) .LT. AMIN) THEN 
AM1N=0UT(INT((N+1)/2+255),I) 

ENDIF 

40 CONTINUE 

DO 50 N-515,1024,2 

0UT(INT((N-:)/2-256),I)=ABS(FT(N))**?/(2*M+l) 
IF (0UT(INT((N-l)/2-256),I) • GT. AMAX) THEN 
AMAX=OUT(1NT(CN-l)/2-256),1) 

ENDIF 

IF (OUT(INT((N l)/2-256),I) .LT. AMIN) THEN 
AMIN=OUT': INT( CN-1)/2-256) , I) 

ENDIF 

50 CONTINUE 

ELSE 

DO 51 N=l,512 

OUT(N,I)=ABS(FT(N))**!/(2*M+1) 

IF (OUT(N,I) .GT. AMAX) THEN 
AMAX=OUTCN,I) 

end;; 
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IF (OUT(N,I) .LT. AMIN) THEN 
AMIN=OUT(N,I) 

ENDIF 

51 CONTINUE 

ENDIF 

10 CONTINUE 

IF (CONTR .EQ. 'C') THEN 
GO TO 98 
ENDIF 

CALL AREA2D(8. ,9. ) 

CALL V0LM3D(10. ,10. ,8. ) 

CALL HEADINCTTL,100,1.,3) 

CALL HEADIN(SIGNAL,100,1.,3) 

CALL HEADINCNO TIME SMOOTHING, 1024 FFT$‘,100,1. ,3) 
CALL MESSAG(WTYPE,100,2.5,9.3) 

CALL INTNO(2*M+1,'ABUT','ABUT') 

CALL MESSAGC P0INTS)$',100,'abut','abut') 

CALL X3NAME('FREQUENCY AXIS$',100) 

CALL Y3NAMF.('TIME AXIS$',100) 

CALL Z3NAME(' $',100) 

CALL VUANGLC-65. ,70. ,700. ) 

CALL XNONUM 
C CALL ZNONUM 

CALL MXlALf ( ' STANDARD' ,'//') 

CALL MX2ALF('L/CGREEK','+') 

CALL ANGLE(-25.0) 

IF ( MODE .EQ. 1 ) THEN 
CALL MESSAGC +0# ',5,0. ,2.3) 

ELSE 

CALL MESSAGC +-P# ',6,0. ,2. 3) 

ENDIF 

CALL ANGLE(-25.0) 

CALL MESSAGC +P# ',5,4.9,0.15) 

CALL GRAF3D(-256,256,256,0,L,L,1. 0*AMIN, 

+ 0.5*(AMAX+AMIN),1.0*AMAX) 

CALL SURMAT(OUT,512,512,1,L,0. ) 

CALL ENDPL(O) 

C 

98 IF (CONTR .NE. 'A') THEN 

CALL CONTOR(OUT,512,L,TTL,SIGNAL,WTYPE,Z,AMAX) 


ENDIF 

111 

CONTINUE 


CALL DONEPL 

99 

STOP 


END 


SSSSSSSSSSSFSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE CONTOR(A,NX,NY,TITLE,SIGNL,VNDW,WLEN,AMAX) 

THIS SUBROUTINE CONTOURS AN NX BY NY ARRAY OF REGULARLY SPACED POINTS. 
NOTE: THE ARRAY MUST BE REAL*4. 

A : SINGLE PRECISION NX BY NY ARRAY OF REGULARLY SPACED POINTS 
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NX: NUMBER OF POINTS IN THE X-DIRECTION 
NY: NUMBER OF POINTS IN THE Y-DIRECTION 
ZINC: CONTOUR INTERVAL 

DIMENSION A(NX,NY) 

COMMON WORKC50000) 

SET PARAMETERS FOR AXES: 

X0RIG--256. 

XSTP=256. 

XMAX=256. 

YORIG=0. 

YSTP=NY 

YMAX=NY 

SET CONTOUR LEVEL 
ZINC=AMAX/10. 

CALL SETCLR('CYAN’) 

SET PAGE AND PLOT SIZES, SET UP AXES FOR PLOT: 

CALL PAGE(8.5,11,0) 

CALL BCOMONC50000) 

CALL AREA2D(6.0,8. 0) 

LABEL AXES: 

CALL XNAMECFREQUENCY - AXIS $',100) 

CALL YNAMECTIME - AXIS $’,100) 

CALL GRAF(XORIG,XSTP,XMAX,YORIG,YSTP,YMAX) 

CALL FRAME 

TITLE: 

CALL HEADINCCONTOUR PLOT$',100,1.,3) 

CALL HEADINCTITLE,100,1. ,3) 

CALL HEADINCSIGNL,100,1. ,3) 

CALL ANGLECO. 0) 

CALL MESSAG(WNDW,100,1. 5,-.7) 

CALL INTNOCWLEH ,’ABUT’,'ABUT') 

CALL MESSAGC P0INTS)$’,100,'abut',’abut') 

MAKE CONTOURS AND DRAW: 

CALL SETCLR('RED’) 

CALL CONMIN(3. 0) 

CALL CONANG(60.) 

CALL CONLINC 0,’MYCON’,’NOLABELS’,2,10) 

CALL CONMAK(A,NX,NY,ZINC) 

CALL CONTURC1,’LABELS’,’DRAW’) 

CALL ENDPL(O) 

RETURN 

END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 
SUBROUTINE MYCON(RARAY,IARAY) 
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DIMENSION RARAY(2),IARAY(1) 

THIS ROUTINE MAKES NEGATIVE CONTOURS DASHED AND THE ZERO LINE HEAVIER. 

CALL RESET('DASH') 

IF (RARAY(l) .GE. 0.) GO TO 10 
CALL DASH 
10 RARAY(2) = 1. 
lARAY(l) = 1 

IF (RARAY(l) .EQ, 0.) lARAY(l) = 2 

RETURN 

END 

ISSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 


* 


** 

* 

* 

CALL FFT(N,XTEMP,X) 

* 

* 

X - OUTPUT COMPLEX ARRAY CONTAINING FFT (1024) 

* 

* 

N - NUMBER OF POINTS 

* 

* 

XTMP - COMPLEX ARRAY CONTAINING DATA SAMPLES 

* 


(starting at l,up to 1024) 

* 

*>V*^V** VoV*iiV*7V * Vc * a Vc***iV Vf WiVVc Vf Vf Vf 



SUBROUTINE FFT(N,XTMP,X) 


FFT00130 


COMPLEX X( 1024) ,XTMP( 1024) ,WTFAC,rt'lP 


FFT00140 


M=INT(L0G10(FLOAT(N))/L0G10(2. )+0. 5) 




EN = N 


FFT00210 


PI = 4. 0*ATAN(1. 0) 


FFT00270 


DO 10 K=0,N-1 


FFT00320 


NEWADR = 0 


FFT00330 


MADDR = K 


FFT00340 


DO 20 1=0,M-1 


FFT00350 


LRMNDR = MOD(MADDR,2) 


FFT00360 


NEWADR = NEWADR + LRMNDR*2**(M-1-I) 


FFT00370 


MADDR = MADDR/2 


FFT00380 

20 

CONTINUE 


FFT00390 


X(NEWADR+1) = XTMP(K+1) 


FFT00400 

10 

CONTINUE 


FFT00410 


DO 50 L=1,M 


FFT00530 


ISPACE = 2**L 


FFT00610 


S = N/ISPACE 


FFT00620 


IWIDTH = ISPACE/2 


FFT00630 


DO 40 J=0,(IWIDTH-1) 


FFT00670 


R = S*J 


FFT00720 


ALPHA = 2.*PI*R/EN 


FFT00730 


WTFAC = CMPLX( COS(ALPHA), -SIN(ALPHA)) 


FFT00740 


DO 30 ITOP=J,N-2,ISPACE 


FFT00750 


IBOT = ITOP + IWIDTH 


FFTOOSGO 


TMP = X(IBOT+l)*WTFAC 


FFT00810 


X(IBOT+l) = X(ITOP+l) - TMP 


FFT00820 


X(ITOP+l) = X(ITOP+l) + TMP 


FFT00830 

30 

CONTINUE 


FFT00840 

40 

CONTINUE 


FFT00850 

50 

COOTINUE 


FfT00860 
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RETURN 

END 


FFTOIOOO 

FFTOlOlO 
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3, 


IPS 


C** THIS FORTRAN FILE COMPUTES THE IPS OF A DATA SEQUENCE *** 

* 

* 

INPUT DATA SEQUENCE IS READ USING FILEDEF 4, AS THE COMPLEX * 
ARRAY X(N). * 

■Z' 

N IS THE LENGTH OF THE DATA SEQUENCE AND IS ADJUSTED FROM THE * 

PARAMETER STATEMENT. N MUST NOT EXCEED 128. * 

* 

ANALYSIS PARAMETERS ARE READ USING FILEDEF 41. THE PARAMETERS * 

ARE: * 

ARGUMENT TYPE ALLOWED VAULUES * 

* 

MODE “II 1 PLOT 0 TO PI * 

2 PLOT PI TO PI * 

* 

PLTR “II 0 SHERPA LASER PRINTER * 

1 IMB79 GRAPHICS TERMINAL * 

■ * 

BWLEN“ 13 3 DIGIT INITIAL WINDOW LENGTH, * 

MUST BE AN ODD INTEGER * 

EWLEN“ 13 3 DIGIT FINAL WINDOW LENGTH * 

WINC “12 2 DIGIT WINDOW INCREMENT, MUST * 

BE AN EVEN INTEGER * 

* 

WTYPE- A19 19 CHARACl'ER STRING USED IN THE * 

PLOT HEADER DISCRIBING THE * 
WINDOW USED. THE CURRENT * 
WINDOW LENGTH IS AUTOMATICALLY * 
INCLUDED * 

* 

CONTR- A1 1 CHARACTER STRING INDICATING * 

TYPE OF PLOT DESIRED * 

A AMPLITUDE PLOT ONLY * 

C CONTOUR PLOT ONLY * 

B BOTH AMPLITUDE AND CONTOUR * 

* 

TTL - A43 43 CHARACTER STRING USED IN THE * 

HEADING WHICH DESCRIBES THE * 

ALGORITHM AND THE CLASS OF * 

SIGNAL USED * 

if 

SIGNAL- A43 45 CHARACTER STRING DESCRIBING * 

TEST SIGNAL 

it 

TSMTH- A25 25 CHARACTEP^ STRING DESCRIBING 

TIPE OF TIME SMOOTHING USED * 

it 
it 

OUT REAL * 

OUTPUT ARRaV of DitlENSlON .M2 BY N * 

it 

SAMp COMPLEX * 
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SHIFTED VERSION OF X 
SAM COMPLEX 

SHIFTED AND CONJUGATED VERSION OF X 
C COMPLEX 

ARRAY OF SUM OF PRODUCTS OF DIMENSION 1 BY 1025 
POSITIONS 1 - 512 ARE CONJUGATE SYMMETRIC WITH 
POSITIONS 514 - 1025. POSITION 513 = (0.,0.) 

FT COMPLEX 

ARRAY OF THE 1024 POINT TRANSFORM C 

Z INTEGER 

LENGTH OF CURRENT WINDOW 

M INTEGER 

MID-POINT OF THE CURRENT WINDOW 
AMAX REAL 

MAXIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS 
AMIN REAL 

MINIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS 

NX HORIZONTAL DIMENSION OF OUT WHICH IS ALWAYS 512 

DATA CAN BE OUTPUT USING FILEDEF 61. POINTS OF INTEREST MUST 
BE DEFINED IN THE APPROPRAITE SECTION OF CODE. BECAUSE 
OF SPACE CONSTRAINTS, THE DATA OUTPUT FILE IS WRITTEN TO 
THE B DISK. 


* 

if 

* 

if 

* 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 


if 

if 


kififififififififififififififififififififififififif‘ 


:itififififififififififififififififififififififififififififififif 


PARAMETER(N= 64) 

COMPLEX X(N),C(1025),SAMP,SAM,FT(1024) 

REAL 0UT(512,N),AMAX,AMIN 

INTEGER NX,K,I,J,MODE,Z,M,BWLEN,EWLEN,WINC,PLTR 
CHARACTER WTYPE*19,TTL*43,SIGNAL*43,TSMTH*25,CONTR*l 
CALL EXCMSC'FILEDEF 4 DISK TEST IN (PERM') 

CALL EXCMSC'FILEDEF 41 DISK PARAM IN (PERM') 

CALL EXCMSC'FILEDEF 61 DISK DATA OUT B (PERM') 

. READ IN PARAMETER LIST . 

READ(41,400)MODE,PLTR,BWLEN,EWLEN,WINC,WTYPE,CONTR,TTL, 

+ STGNAL,TSMTH 

400 FORMAT (1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X,A43 
+ /:X,A43/1X,A25) 


.... test to ENSURE WINDOW LENGTH IS APPROPRIATE 
I=N-1 

IF ((BWLEN .GT. I) .OR. (EWLEN .GT. I)) THEN 

f TT^ Y 'rr* / \ 

yUy j 

GO TO 99 
ENDIF 
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) 


I=M0D(BWLEN,2) 

K=M0D(WINC,2) 

IF (I .EQ. 0) THEN 

IF (K .EQ. 1) THEN 
WRITE(*,68) 

GO TO 99 

ELSE 

WRITE(*,67) 

GO TO 99 
ENDIF 

ENDIF 

69 FORMAT (IX,'WINDOW LENGTO EXCEEDS LENGTH OF THE DATA') 

68 FORMAT (IX,'WINDOW INCREMENT MUST BE EVEN') 

67 FORMAT (IX,'INITIAL WINDOW LENGTH MUST BE ODD') 


-—PLOTTING DEVICE CALL 
IF (PLTR.EQ.0)THEN 
CALL COMPRS 

ELSE 

CALL IBM79 
ENDIF 


PI=4*ATAN(1. ) 

NX=512 

READ (4,*)(X(J),J=1,N) 

DO 111 Z= BWLEN,EWLEN,WINC 
WRITE(61,600)TTL,SIGNAL,TSMTH,WTYPE,Z 
600 FORMAT (1X,A43/1X,A43/1X,A25/1X,A19,I3,' POINTS)') 

CALL ANGLE(0. ,0. ) 

M=(Z-l)/2 

ANAX=0. 

AMIN=4MAX 
DO 10 1=1.N 
DO 20 K=0,512 
SAMP=(0.,0.) 

SAM=SAMP 

IF ( (I+K) .LE. N ) THEN 
SAMP=X(I+K) 

ENDIF 

IF ( (I-K) .GT. 0 ) THEN 
SAM=CONJG(X(I-K)) 

ENDIF 

C 

C+ - + - + - + - + - + - + - + - + - + - + - + - + - + - + “ + “ + - + 
C+-+-+-+ sum of product is IPS difference of product is IMRD -+-+-+ 
IF (K .LE. M) THEN 

C(K+1)=(X(I)*SAM + CONJG(X(I))*SAMP) 

+ *(0. 54+0. 46*CO.sr2*PI*K/(2*M))) 

ELSE 

C(K+1)=0 

hiNUiJr 

C+ -+-+-+-+-+-+-+-+-+-+-+-+- +-+-+-+-+ 
C 
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C(1024-K+1)=CONJG(C(K+1)) 

20 CONTINUE 

C(513)=(0. ,0. ) 

CALL FFT(1024, :,FT) 

DO 40 K=1,M0DE*j12,M0DE 
IF (REAL(FT(K)) . GT. AMAX) THEN 
AMAX=REAL(FT(K)) 

ENDIF 

IF (REAL(FT(K)) .LT. AMIN) THEN 
AMIN=REAL(FT(K)) 

ENDIF 

IF ( (K .LT; 514) .AND. (MODE . EQ. 2)) THEN 
OUT(INT((K+1)/2+255),I)=REAL(FT(K)) 

ELSE 

IF ( MODE .EQ. 2 ) THEN 

OUT(INT((K+l)/2-257),I)=REAL(FT(K)) 

ELSE 

OUT(K,I)=REAL(FT(K)) 

ENDIF 

ENDIF 

40 CONTINUE 

10 CONTINUE 

C.FOR TIME SMOOTHING PURPOSES . 

DO 48 K=l,512 
DO k( 1=1,N-2 

OUr(K,I)=(OUT(K,I)+OUT(K,I+l)+OUT(K I+2))/3 

46 CONTINUE 

DO 47 I=N,3,-1 

OUT(K,I)=(OUT(K,I)+OUT(K,I-l)+OUT(K,I-2))/3 

47 CONTINUE 

. DATA OUTPUT . 

IF((K.GE. 208).AND.(K. LT.370))THEN 
IF(((K.GE.120).AND.(K.LT.140)).OR.((K.GE.370).AND. 
+ (K, LT. 390)))THEN 

WRITE(61,601) 

DO 81 I=1,N-1,14 

WRITE(61,602)K,I,0UT(K,I) 

81 CONTINUE 

ENDIF 

601 FORMAT ('FREQ BIN=',8X,'TIME BIN=',7X,'AMPLITUDE=') 

602 FORMAT (10X,I4,13X,I4,14X,E14. 7) 


48 


CONTINUE 


. PLOTTING 

IF(CONTR.EQ.'C')THEN 
GO TO 50 
ENDIF 

CALL BSHIFT ( -0.2 ,-.25) 
CALL AREA2D(8,9) 

CALL VOLM3D(10,10,8) 

CALL HEADINCTTL,100,1. ,3) 
CALL HEADINCSIGNAL, 100,1., ,3) 
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CALL HEADIN(TSMTH,100,1. ,3) 

CALL MESSAG( W'fYPE, 100,2. 5,9, 3) 

CALL INTiMO(Z,'ABUT' 'ABUT*) 

CALLMESSAGC' POINT.. )$',100 ’ABUT’,* ABUT') 

CALL X3NAME('FREQUENCY AXIS$',100) 

CALL y3NAME(’TIME AXIS$',106) 

CALL Z3NAME(' $*,100) 

CALL VUANGL(-65,70,700) 

CALL XNONUM 
C CALL ZNONUM 

CALL MX1ALF('STANDARD','#') 

CALL MX2ALF(' L/CGREEK’ 

CALL ANGLE(-25.0) 

IF ( MODE .EQ. 2 ) THEN 

CALLMESSAGC +-P# ',6,0. ,2.3) 

ELSE 

CALL MESSAGC' +0# ',5,0.,2.3) 

ENDIF 

CALL ANGLE(-25. 0) 

CALLMESSAGC +P# ’,5,4.9,0.15) 

CALL GRAF3D(-256,256,256,1,N,N,1. 0*AMIN, 

+ 0.5*(AMAX-AMIN),1. 0*AMAX) 

CALL SURHAT(0UT,512,512,1,N,0.) 

CALL ENDPL(O) 

C 

50 IF (CONTR.NE.'A')THEN 

DO 49 1=1,N 

DO 51 K=l,512 

IF (0UT(K,I) .LT. 0) THEN 
0UT(K,I)=0 
ENDIF 

51 CONTINUE 

49 CONTINUE 

CALL CONTOR(OUT,NX,N,TTL,SIGNAL,WTYPE,TSMTH,Z,AMAX) 
ENDIF 


WRITE(61,603)AMAX,AMIN 
603 FORMAT (IX,'MAXIMUM AMPLITUDE=',E14.7 
+ /'MINIMUM AMPLITUDE=',E14. 7) 

111 CONTINUE 

CALL DONEPL 
99 STOP 
END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE CONTOR(A,NX,NY,TITLE,SIGNL,WNDW,TAVG,WLEN,AMAX) 

THIS SUBROUTINE CONTOURS AN NX BY NY ARRAY OF REGULARbY SPACED POINTS. 
NOiE: THE ARRAY MUST BE REAL*4. 

A : SINGLE PRECISION NX BY NY ARRAY OF REGULARLY SPACED POINTS 
NX: NUMBER OF POINTS IN THE X-DIRECTION 
NY: NUMBER OF POINTS IN THE Y-DIRECTION 
ZINC: CONTOUR INTERVAL 
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DIMENSION A(NX,NY) 

COMMON WORK(50000) 

SET PARAMETERS FOR AXES: 

X0RIG=-256. 

XSTP=256. 

XMAX=256. 

Y0RIG=0. 

YSTP=NY 

YMAX=NY 

SET CONTOUR LEVEL 
ZINC=AMAX/10. 

CALL SETCLR(’CYAN') 

SET PAGE AND PLOT SIZES, SET UP AXES FOR PLOT: 

CALL PAGE(8. 5,11.0) 

CALL BCOMONC50000) 

CALL AREA2D(6.0,8. 0) 

LABEL AXES: 

CALL XmMECFREQUENCY - AXIS $',100) 

CALL YNAMECTIME - AXIS $’,100) 

CALL GRAF(XORIG,XSTP,XMAX,YORIG,YSTP,YMAX) 

CALL FRAME 

TITLE: 

CALL HEADINCCONTOUR PL0T$',100,1. ,4) 

CALL HEADINCTITLE,100,1. ,4) 

CALL HEADINCSIGNL,100,1. .4) 

CALL HEADINCTAVG,100,1. ,4) 

CALL ANGLE(O.O) 

CALL MESSAG(WNDW,100,1.5,-.7) 

CALL INTNOCWLEN ,'ABUT','ABUT') 

CALL MESSAGC POINTS)?',100,'ABUT'abut') 

MAKE CONTOURS AND DRAW: 

CALL SETCLR('RED') 

CALL CONMINO, 0) 

CALL C0NANG(60.) 

CALL CONLIN(0,'MYCON','NOLABELS',2,10) 

CALL CONMAK(A,NX,NY,ZINC') 

CALL CONTURC1,'LABELS','DRAW’) 

CALL ENDPL(O) 

RETURN 

END 

>SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE MYC0N(RARAY,IARAY) 

DIMENSION RARAY(2),IARAY(1) 

THIS ROUTINE MAK. NEGATIVE CONTOURS DASHED AND THE ZERO LINE HEAVIER. 
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CALL RESET('DASH’) 

IF (RARAY(l) .GE. 0.) GO TO 10 
CALL DASH 
10 RARAY(2) = 1. 
lARAY(l) = 1 

IF (RARAY(l) .EQ. 0.) lARAY(l) = 2 
RETURN 
END 
C 

csssssssssssssssssssssssssssssssssss 


it * 

* CALL FFT(N,XTEMP,X) * 

* * 

* X - OUTPUT COMPLEX ARRAY CONTAINING FFT (1024) * 

* N - NUMBER OF POINTS * 

* XTMP - COMPLEX ARRAY CONTAINING DATA SAMPLES * 

* (starting at ljUp to 1024) * 





SUBROUTINE FFT(N,XTMP,X) 

FFT00130 


COMPLEX X(1024),XTMP(1024),WTFAC,TMP 
M=INT(LOG10(FLOAT(N))/LOG10(2.)+0.5) 

FFT00140 


EN = N 

FFT00210 


PI = 4.0*ATAN(1. 0) 

FFT0027 


DO 10 K=0,N-1 

FFT00320 


NEWADR = 0 

FFT00330 


MADDR = K 

FFT00340 


DO 20 1=0,M-1 

FFT00350 


LRMNDR = MOD(MADDR,2) 

FFT00360 


NEWADR = NEWADR + LRMNDR*2**(M-1-I) 

FFT00370 


MADDR = MADDR/2 

FFT00380 

20 

CONTINUE 

FFT00390 


X(NEWADR+1) = XTMP(K+1) 

FFT00400 

10 

CONTINUE 

FFT00410 


DO 50 L=1,M 

FFT00530 


ISPACE = 2**L 

FFT00610 


S = N/ISPACE 

FFT00620 


IWIDTH = ISPACE/2 

FFT00630 


DO 40 J=0,(IWIDTH-1) 

FFT00670 


R = S*J 

FFT00720 


ALPHA = 2.*PI*R/EN 

FFT00730 


WTFAC = CMPLX( COS(ALPHA), -SIN(ALPHA)) 

FFT00740 


DO 30 ITOP=J,N-2,ISPACE 

FFT00750 


IBOT = ITOP + IWIDTH 

FFT00800 


TMP = X(IBOT+l)*WTFAC 

FFT00810 


X(IBOT+l) = X(ITOP+l) - TMi"- 

FFT00820 


X(ITOP+l) = X(ITOP+l) + TMP 

FFT00830 

30 

CONTINUE 

FFT00840 

40 

CONTINUE 

FFT00850 

50 

CONTINUE 

FFT00860 


RETURN 

FFTOIOOO 


END 

FFTOlOlO 
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4. PWD 


c** this fortran file computes the PWD OF A DATA SEQUENCE *** 

* 
* 

INPUT DATA SEQUENCE IS READ USING FILEDEF 4, AS THE COMPLEX * 
ARRAY X(N). * 

* 

N IS THE LENGTH OF THE DATA SEQUENCE AND IS ADJUSTED FROM THE * 
PARAMETER STATEMENT. N MUST NOT EXCEED 128. * 

* 

ANALYSIS PARAMETERS ARE READ USING FILEDEF 41. THE PARAMETERS * 
ARE: * 

ARGUMENT TYPE ALLOWED VAULUES * 

* 

MODE - II 1 PLOT 0 TO PI * 

2 PLOT PI TO PI * 

it 

PLTR - II 0 SHERPA LASER PRINTER * 

1 IMB79 GRAPHICS TERMINAL * 

■ * 

BWLEN- 13 3 DIGIT INITIAL WINDOW LENGTH, * 

MUST BE AN ODD INTEGER * 

EWLEN- 13 3 DIGIT FINAL WINDOW LENGTH * 

WINC - 12 2 DIGIT WINDOW INCREMENT, MUST * 

BE AN EVEN INTEGER 

it 

WTYPE- A19 19 CHARACTER STRING USED IN THE * 

PLOT HEADER DISCRIBING THE * 

WINDOW USED. THE CURRENT * 

WINDOW LENGTH IS AUTOMATICALLY * 
INCLUDED * 

* 

CONTR- A1 1 CHARACTER STRING INDICATING * 

TIPE OF PLOT DESIRED * 

it 

A AMPLITUDE PLOT ONLY * 

C CONTOUR PLOT ONLY * 

B BOTH AMPLITUDE AND CONTOUR * 

* 

TTL - A43 43 CHARACTER STRING USED IN THE * 

HEADING WHICH DESCRIBES THE * 

ALGORITHM AND THE CLASS OF * 

SIGNAL USED * 

it 

SIGNAL- A43 43 CHARACTER STRING DESCRIBING * 

TEST SIGNAL * 

it 

TSMTH- A25 25 CHARACTER STRING DESCRIBING * 

TYPE OF TIME SMOOTHING USED * 

it 
it 

OUT REAL * 

OUTPUT ARRAY OF DIMENSION 512 BY N 

it 

SAMP COMPLEX * 
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SHIFTED VERSION OF X * 

* 

SAM COMPLEX * 

SHIFTED AND CONJUGATED VERSION OF X * 

* 

C COMPLEX * 

ARRAY OF SUM OF PRODUCTS OF DIMENSION 1 BY 1025 * 

POSITIONS 1-512 ARE CONJUGATE SYMMETRIC WITH * 

POSITIONS 514 - 1025. POSITION 513 = (0.,0.) * 

* 

FT COMPLEX * 

ARRAY OF THE 1024 POINT TRANSFORM C * 

* 

Z INTEGER * 

LENGTH OF CURRENT WINDOW * 

* 

M INTEGER * 

MID-POINT OF THE CURRENT WINDOW * 

* 

AMAX REAL * 

MAXIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS * 

* 

AMIN REAL * 

MINIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS * 

* 

NX HORIZONTAL DIMENSION OF OUT WHICH IS ALWAYS 512 * 

* 


DATA CAN BE OUTPUT USING FILEDEF 61. POINTS OF INTEREST MUST * 
BE DEFINED IN THE APPROPRAITE SECTION OF CODE. BECAUSE * 

OF SPACE CONSTRAINTS, THE DATA OUTPUT FILE IS WRITTEN TO * 

THE B DISK. * 




PARAMETER(N= 64) 

COMPLEX X(512),C(1025),SAMP,SAM,FT(1024) 

REAL OUT(512,N),AMAX,AMIN 

INTEGER NX,K,I,J,MODE,Z,M,BWLEN,EWLEN,WINC,PLIR 
CHARACTER WTYPE*19,TTL*43,SIGNAL*43,TSMTK*25.C0NTR*1 
CALL EXCMS('FILEDEF 4 DISK TEST IN (PERM') 

CALL EXCMS('FILEDEF 41 DISK PARAM IN (PERM') 

CALL EXCMS('FILEDEF 61 DISK DATA OUT B (PERM') 

.. READ IN THE PARAMETER LIST. 

READ(41,400)MODE,PLTR,BWLEN,EWLEN,WINC,WTYPE,CONTR,TTL, 

+ SIGNAL,TSMTH 

400 F0RMAT(1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X,A43 
+ /1X,A43/1X,A25) 


. TEST TO ENSURE WINDOW LENGTH IS APPROPRIATE 

1=N-1 

IF ((BWLEN .GT. I) .OR. (EWLEN .GT. I)) THEN 
WRITE(*,69) 

GO to'99 
ENDIF 
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I«M0D(BWLEN.2'» 

K-M0D(WINC,2) 

IF (I .EQ. 0) THEF 

IF (K . EQ. 1) THEN 
WRITE(*,68) 

GO TO 99 

ELSE 

WRrfE(*,67) 

GO TO 99 
ENDIF 

ENDIF 

69 FORMAT fIX,'WINDOW LENGTH EXCEEDS LENGTH OF THE DATA') 

68 FORMAT (IX,'WINDOW INCREMENT MUST BE EVEN') 

67 FORMAT (IX,'INITIAL WINDCW LENGTH MUST BE ODD') 


■.PLOTTING DEVICE CALL 

IF(PLTR.EQ. 0)THEN 
CALL COMPRS 
ELSE 

CALL IBM79 
ENDIF 


PI=4*ATAN(1. ) 

NX=512 

READ(4,*)(X(J),J=1,2*N,2) 

. ..DATA INTERPOLATION 

DO 5 J=2,2*N,2 
X(J)=(0, ,0.) 

5 CONTINUE 

CALL FFT(2*N,X,FT) 

DO 10 J=N/2+2,2*N-N/2+l 
FT(J)=(0. ,0.) 

10 CONTINUE 

DO 20 J=1,2*N 

FT(J)=CONJG(FT(J)) 

20 CONTINUE 

CALL FFT(2*N,FT,X) 

DO 30 J=1,2*N 

X(J)=CONJG(X(J))/N 
30 CONTINUE 


DO 111 Z=BWLEN,EWLEN,WINC 
WRITE(61,600)TTL,SIGNAL,TSMTH,WTYPE,Z 
600 F0RMAT(1X,A43/1X,A37/1X,A25/1X,A19,I3,’ POINTS)') 

M=(Z-l)/2 
AMAX=0 

DO 40 I=1,2*N,2 
DO 50 K=0,512 
SAMP=(0. ,0. ) 
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SAM=(0.,0.) 

IF ( (I+K) .LE. 2*N ) THEN 
SAMP=X(I+K) 

ENDIF 

IF ( (I-K) .GT. 0 ) THEN 
SAM=CONJG(X(I“K)) 

ENDIF 

C 

C 

IF (K .LE. 2*M) THEN 

C(K+1)=SAMP*SAM*(0.54+0.46*C0S(2*PI*K/(4*M))) 
ELSE 

C(K+1)=0 

ENDIF 

C 

C 

C(1024-K+1)=C0NJG(C(K+1)) 

50 CONTINUE 

C(513)=(0. ,0. ) 

CALL FFT(1024,C,FT) 

DO 60 K=1,M0DE*512,M0DE 

IF (REAL(FT(K)) . GT. AMAX) THEN 
AMAX=REAL(FT(K)) 

ENDIF 

IF (REAL(FT(K)) . LT. AMIN) THEN 
AMIN=REAL(FT(K)) 

ENDIF 

IF ( (K.LT.514) .AND. (MODE .EQ. 2)) THEN 
OUT(INT((K+1)/2+255),(I+l)/2)=REAL(FT(K)) 

ELSE 

IF( MODE .EQ. 2 ) THEN 
OUT(INT((K+1)/2-257),(I+l)/2)=REAL(FT(K)) 
ELSE 

0UT(K,(I+1)/2)=REAL(FT(K)) 

ENDIF 
ENDIF 
60 CONTINUE 

40 CONTINUE 
C 

C.FOR TIME SMOOTHING PURPOSES. 

DO 48 K=l,512 
DO 46 1=1,N-2 

OUT(K,I)=(OUT(K,I)+OUT(K,I+l)+OUT(K,I+2))/3 

46 CONTINUE 

DO 47 I=N,3,-1 

0UT(K,I)=(0UT(K,I)+0UT(K,I-l)+0UT(K,I-2))/3 

47 CONTINUE 
C 

C.DATA OUTPUT. 

IF((K. GT. 285). AND. (K. LT. 310))THEN 
WRITE(61,601) 

DO 81 1=1,N 

WRITE(6l,602)K,I,OUT(K,I) 

81 CONTINUE 
ENDIF 

601 FORMATCFREQ BIN=',8X,'TIME BIN=',7X,'AMPLITIDE=') 
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602 FORMAT(10X,I4,13X,I4,14X,E14.7) 
C. 

c 

48 CONTINUE 


■. PLOTTING . 

IF (CONTR.EQ.'C')THEN 
GO TO 52 
ENDIF 

CALL BSHIFT ( -0.2 , -0.25) 

CALL AREA2D(8,9) 

CALL VOLM3D(10,10,8) 

CALL HEADIN(TTL,100,1.,3) 

CALL HEADINCSIGNAL,100,1. ,3) 

CALL HEADINCTSMTH,100,1.,3) 

CALL ANGLE(0. ,0. ) 

CALL MESSAG(WTYPE,100,2.5,9.3) 

CALL INTNO(Z ,'ABUT*,'ABUT') 

CALL MESSAGC points )$*,100,'ABUT*,'ABUT') 

CALL X3NAME('FREQUENCY AXIS$S100) 

CALL Y3NAME('TIME AXIS$',100) 

CALL Z3NAME(' $’,100) 

CALL VUANGLC-65,70;700) 

CALL XNONUM 
C CALL ZNONUM 

CALL MX1ALF('STANDARD','#') 

CALL MX2ALF('L/CGREEK*,'+') 

CALL ANGLE(-25.0) 

IF ( MODE .EQ. 2 ) THEN 
CALL MESSAGC +-P// *,6,0. ,2. 3) 

ELSE 

CALL MESSAGC* +0# ',5,0.,2.3) 

ENDIF 

CALL ANGLE(-25.0) 

CALL MESSAGC* +PC/ ',5,4.9,0.15) 

CALL GRAF3DC-256,256,256,1,N,N,1.0*AMIN,. 5*CAMAX-AMIN), 

+1.0*AMAX) 

CALL SURMATCOUT,512,512,1,N,0.) 

CALL ENDPLCO) 

C 

52 IFCCONTR. NE. 'A*)THEN 
DO 49 1=1,N 

DO 51 K=l,512 

IFCCOUTCK,I).LT. 5.0).AND.C0UTCK,I) . GT. -5.0))THEN 
0UTCK,I)=0 
ENDIF 

51 CONTINUE 

49 CONTINUE 

CALL CONTORC OUT,NX,N,TTL,SIGNAL,WTYPE,TSMTH,Z) 

ENDIF 

C. 

c 

WRITEC61,603)AMAX,AMIN 
603 FORMATCIX,'MAXIMUM AMPLITUDE=',E14. 7 
+ /IX,'MINIMUM AMPLITUDE=',E14. 7) 
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111 


CONTINUE 
CALL DONEPL 
99 STOP 
END 

ssssssssssssssssssssssssssssssssss 

SUBROUTINE CONTOR(A,NX,NY.TITLE,SIGNL,WNDW,TAVG,WLEN) 

THIS SUBROUTINE CONTOURS AN NX BY NY ARRAY OF REGULARLY SPACED POINTS. 
NOTE: THE ARRAY MUST BE REAL*4. 

A : SINGLE PRECISION NX BY NY ARRAY OF REGULARLY SPACED POINTS 
NX: NUMBER OF POINTS IN THE X-DIRECTION 
NY; NUMBER OF POINTS IN THE Y-DIRECTION 
ZINC: CONTOUR INTERVAL 

DIMENSION A(NX,NY) 

COMMON WORK(50000) 

SET PARAMETERS FOR AXES: 

XORIG=-256. 

XSTP=256. 

XMAX=256. 

YORIG=0. 

YSTP=NY 

YMAX=NY 

SET CONTOUR LEVEL 
ZINC=AMAX/10 

CALL SETCLR('CYAN') 

SET PAGE AND PLOT SIZES, SET UP AXES FOR PLOT: 

CALL PAGE(8.5,11.0) 

CALL BCOMONC50000) 

CALL AREA2D(6.0,8. 0) 

LABEL AXES: 

CALL XNAME('FREQUENCY - AXIS $',100) 

CALL YNAMECTIME - AXIS $’,100) 

CALL XNONUM 

CALL GRAF(XORIG,XSTP,XMAX,YORIG,YSTP,YMAX) 

CALL FRAME 

TITLE: 

CALL HEADINCCONTOUR PLOT$',100,1. ,4) 

CALL HEADINCTITLE,100,1. ,4) 

CALL HEADINCSIGNL,100,1. ,4) 

CALL HEADINCTAVG,100,1.,4) 

CALL ANGLECO.O) 

CALL MESSAGCWNDW,100,1.5.*.7) 

k T T / t tT t N 

Lii41(JLi , aDVi 5 HDVl ) 

CALL MESSAGC'POINTS)$',100,'ABUT','ABUT') 

C 

C MAKE CONTOURS AND DRAW: 
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CALL SETCLR('RED') 

CALL CONMINO.O) 

CALL CONANG(60, ) 

CALL C0NLIN(0,’MYG0N'NOLABELS',2,10) 

CALL CONMAKCA,NX,NY,'SCALE') 

CALL CONTURC1,’LABELS','DRAW') 

CALL ENDPL(O) 

RETURN 

END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE KYCON(RARAY,lARAY) 

DIMENSION R/1RAY(2), IARAY( 1) 

THIS ROUTINE MAKES NEGATIVE CONTOURS DASHED AND THE ZERO LINE HEAVIER. 

CALL RESET('DASH') 

IF (RARAYd) .GE. 0.) GO TO 10 
CALL DASH 
10 RARAY(2) = 1. 
lARAY(l) = 1 

IF (RARAY(l) .EQ. C.) lARAY(l) = 2 

RETURN 

END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 


* 

* 

* 

ic 

Vf 

* 

* 


CALL FFT(N,XTEMP,X) 

X - OUTPUT COMPLEX ARRAY CONTAINING FFT (1024) 
N • NUMBER OF POINTS 

XTMP • COMPLEX ARRAY CONTAINING DATA SAMPLES 
(starting at l,up to 1024) 


it 

* 

* 

* 

* 

* 

* 


SUBROUTINE FFT(N,XTMP,X) 

COMPLEX X(1024),XTMP(1024),WTFAC,TMP 
M=INT(LOG10(FLOAT(N))/LOG10(2.)+0.5) 

EN = N 

PI = 4.C*ATAN(1. 0) 

DO 10 K=0,N-1 
NEWADR = 0 
MADDR = K 
DO 20 1=0,M-1 

LRMNDR = M0D(MADDR,2) 

NEWADR = NEWADR + LRMNDR*2**(M-1-I) 
MADDR = MADDR/2 
20 CONTINUE 

X(NEWADR+1) = XTMP(K+1) 

10 CONTINUE • ■ 

DO 50 L=1,M 

ISPACE = 2**L 


FFT00130 

FFT00140 

FFT00210 

FFT00270 

FFT00320 

FFT00330 

FFT00340 

FFT00350 

FFT00360 

FFT00370 

FFT00380 

FFT00390 

FFT00400 

FFT00410 

FFT00530 

FFT00610 
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S = N/ISPACE 

FFT00620 


IWIDTH = ISPACE/2 

FFT00630 


DO 40 J-0,(IWIDTH-1) 

FFT00670 


R = S*J 

FFT00720 


ALPHA = 2.*PI*R/EN 

FFT00730 


WTFAC = CMPLXC COS(ALPHA), -SIN(ALPHA)) 

FFT00740 


DO 30 IT0P=J,N-2,ISPACE 

FFT00750 


IBOT = ITOP + IWIDTH 

FFT00800 


TMP = X(IB0T+1)*WTFAC 

FFT00810 


X(IB0T+1) = X(IT0P+1) - TMP 

FFT00820 


X(IT0P+1) = X(IT0P+1) + TMP 

FFT00830 

30 

CONTINUE 

FFT00840 

40 

CONTINUE 

FFT00850 

50 

CONTINUE 

FFT00860 


RETURN 

FFTOIOOO 


END 

FFTOlOlO 
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5 . \RD,\^ 


CM, this fortran file computes the squared magnitude of the *** 

C** RIHACZEK DISTRIBUTION *** 

* 

INPUT DATA SEQUENCE IS READ USING FILEDEF 4, AS THE COMPLEX * 

ARRAY X(N). * 

* 

N IS THE LENGT.. OF THE DATA SEQUENCE AND IS ADJUSTED FROM THE * 
PARAMETER STATEMENT. N MUST NOT EXCEED 128. * 

* 

ANALYSIS PARAMETERS ARE READ USING FILEDEF 41. THE PARAMETERS * 
ARE; * 

ARGUMENT TYPE ALLOWED VAULUES * 

* 

MODE - II 1 PLOT 0 TO PI * 

2 PLOT PI TO PI * 

* 

PLTR - II 0 SHERPA LASER PRINTER * 

1 IMB79 GRAPHICS TERMINAL * 

•k 

BWLEN- 13 3 DIGIT INITIAL WINDOW LENGTH, * 

MUST BE AN ODD INTEGER * 

EWLEN- 13 3 DIGIT FINAL WINDOW LENGTH * 

WINC - 12 2 DIGIT WINDOW INCREMENT, MUST * 

BE AN EVEN INTEGER * 

* 

WTYPE- A19 19 CHARACTER STRING USED IN THE * 

PLOT HEADER DISCRIBING THE * 

WINDOW USED. THE CURRENT * 

WINDOW LENGTH IS AUTOMATICALLY * 

INCLUDED * 

* 

CONTR- A1 1 CHARACTER STRING INDICATING * 

TYPE OF PLOT DESIRED 

* 

A AMPLITUDE PLOT ONLY * 

C CONTOUR PLOT ONLY * 


B BOTH AMPLITUDE AND CONTOUR * 

* 

TTL - A43 43 CHARACTER STRING USED IN THE * 

HEADING WHICH DESCRIBES THE * 
ALGORITHM AND THE CLASS OF * 


SIGNAL USED * 

* 

SIGNAL- A43 43 CHARACTER STRING DESCRIBING * 

TEST SIGNAL * 

* 

TSMTH- A25 25 CHARACTER STRING DESCRIBING * 

TYPE OF TIME SMOOTHING USED * 

it 

■It 

OUTR REAL * 

INITIALLY USED AS STORAGE FOR THE TIME SMOOTHED IPS, THEN * 
AS OUTPUT ARRAY OF DIMENSION 512 BY N * 
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* 

OUTI REAL * 

ARRAY USED TO HOLD THE THE TIME SMOOTHED IMRD 

* 

SAMP COMPLEX * 

SHIFTED VERSION OF X * 

* 

SAM COMPLEX * 

SHIFTED AND CONJUGATED VERSION OF X * 

* 

RE COMPLEX * 

ARRAY OF SUM OF PRODUCTS OF DIMENSION 1 BY 1025, * 

POSITIONS 1 - 512 ARE CONJUGATE SYMMETRIC WITH * 

POSITIONS 514 - 1025, POSITION 513 = (0.,0,). * 

* 

IM COMPLEX * 

ARRAY OF DIFFERENCE OF PRODUCTS OF DIMENSION 1 BY 1025, * 

POSITIONS 1 - 512 ARE CONJUGATE SYMMETRIC WITH * 

POSITIONS 514 - 1025, POSITION 513 = (0.,0.)• * 

* 

FTR COMPLEX * 

ARRAY OF THE 1024 POINT TRANSFORM RE * 

* 

FTI COMPLEX * 

ARRAY OF THE 1024 POINT TRANSFORM IM 

•* 

Z INTEGER * 

LENGTH OF CURRENT WINDOW * 

* 

M INTEGER * 

MID-POINT OF THE CURRENT WINDOW * 

* 

AMAX REAL * 

MAXIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS * 

* 

AMIN REAL * 

• MINIMUM AMPLITUDE, USED TO SCALE VERTICAL AXIS * 

* 

NX HORIZONTAL DIMENSION OF OUT WHICH IS ALWAYS 512 * 

* 

DATA CAN BE OUTPUT USING FILEDEF 61. POINTS OF INTEREST MUST * 
BE DEFINED IN THE APPROPRAITE SECTION OF CODE. BECAUSE * 
OF SPACE CONSTRAINTS, THE DATA OUTPUT FILE IS WRITTEN TO * 
THE B DISK. * 

* 


PARAMETER(N= 64) 

COMPLEX X(N),RE(1025),SAMP,SAM,FTR(1024) 

COMPLEX IM(1025),FTI(1024) 

REAL OUTR(512,N),OUTI(512,N),AMAX,AMIN 
INTEGER NX,K,1,J,MODE,Z.M,BWLEN.EWLEN.WING,PLTR 
CHARACTER WTYPE*19,TTL*43,SIGNAL*37,TSM7H*25,CONTR*l 
CALL EXCMS('FILEDEF 4 DISK TEST IN (PERM') 

CALL EXCMS('FILEDEF 41 DISK PARAM IN (PERU ) 
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CALL EXCMSC'FILEDEF 61 DISK DATA OUT B (PERM’) 

■.READ IN PARAMETER LISITNG. 

READ (41,400)MODE,PLTR,BWLEN,EWLEN,WING,WTYPE,CONTR,TTL, 

+ SIGNAL,TSMTH 

FORMAT (1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X,A43 
+ /1X,A37/1X,A2S) 


. TEST TO ENSURE WINDOW LENGTH IS APPROPRIATE . 

I=N-1 

IF ((BWLEN .GT. I) .OR. (EWLEN .GT. I)) THEN 
WRITE(*,69) 

GO TO 99 
ENDIF 

I=MOD(BWLEN,2) 

K=MOD(WINC,2) 

IF (I .EQ. 0) THEN 

IF (K .EQ. 1) THEN 
WRITE(*,68) 

GO TO 99 

ELSE 

WRITE(*,67) 

GO TO 99 
ENDIF 

ENDIF 

69 FORMAT (IX,’WINDOW LENGTH EXCEEDS LENGTH OF THE DATA’) 

68 FORMAT (IX,’WINDOW INCREMENT MUST BE EVEN’) 

67 FORMAT (IX,’INITIAL WINDOW LENGTH MUST BE ODD’) 


. PLOTTING DEVICE CALL 

IF (PLTR.EQ. 0)THEN 
CALL COMPRS 

ELSE 

CALL TBM79 
ENDIF 


PI=4*ATAN(1. ) 

NX=512 

READ (4,*)(X(J),J=1,N) 

DO 111 Z= BWLEN,EWLEN,WINC 
WRITE (61,600)TrL,SIGNAL,TSMTH,WTyPE,Z 
600 FORMAT (1X,A43/1X,A37/1X,A25/1X,A19,I3,’ POINTS)') 
CALL ANGLE(0. ,0. ) 

M=(Z-l)/2 

AMAX=0. 

AMIN=AMAX 
DO 10 1=1,N 
DO 20 K=0,512 
SAMP=(0. ,0. ) 

SAM=SAMP 

IF ( (I+K) .LE. N ) THEN 
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SAMP=X(I+K) 

ENDIF 

IF ( (I-K) .GT. 0 ) THEN 
SAM=CONJG(X(I-K)) 

ENDIF 

C 

IF (K .LE. M) THEN 

RE(K+1)=(X(I)*SAM + CONJG(X(I))*SAMP) 

+ *(0.54+0.46*COS(2*PI*K/(2*M))) 

C 

IM(K+1)=(X(I)*SAM - CONJG(XCI))*SAMP) 

+ *(0.54+0.46*C0S(2*PI*K/(2*M))) 

ELSE 

RE(K+1)=0 

IM(K+1)=0 

ENDIF 

C 

RE(1024-K+1)=CONJG(RE(K+1)) 

IM(1024-K+1)*C0NJG(IM(K+1)) 

20 CONTINUE 

RE(513)=(0.,0.) 

IM(513)=(0.,0,) 

CALL FFT(1024,RE,FTR) 

CALL FFT(1024,IM,FTI) 

DO 40 K=l,MODE*512,MODE 
C 

IF ( (K .LT. 514) .AND. (MODE .EQ. 2)) THEN 
OUTR(INT((K+1)/2+255),I)=REAL(FTR(K)) 

OUTI(INT((K+1)/2+255),I)=REAL(FTI(K)) 

ELSE 

IF ( MODE .EQ. 2 ) THEN 

OUTR(INT((K+1)/2-257),I)-REAL(FTR(K)) 

OUTI(INT((K+1)/2“257),I)=REAL(FTI(K)) 

ELSE 

OUTR(K,I)=REAL(FTR(K)) 

OUTI(K,I)=REAL(FTI(K)) 

ENDIF 

ENDIF 

40 CONTINUE 

10 CONTINUE 

C 

C.FOR TIME SMOOTHING. 

DO 48 K=l,512 
DO 46 1=1,N-2 

OUTR(K,I)=(OUTR(K,I)+OUTR(K,I+l)+OUTR(K,1+2))/3 
OUTI(K,I)=(OUTI(K,I)+OUTI(K,I+l)+OUTI(K,I+2))/3 

46 CONTINUE 

DO 47 I=N,3,-1 

OUTR(K,I)=(OUTR(K,I)+OUTR(K,I-l)+OUTR(K,I-2))/3 

0UTI(K,I)=(0UTI(K,I)+0UTI(K,I-l)+0UTI(K,I-2))/3 

47 CONTINUE 

48 CONTINUE 
C 

C.. . 

C+ - + - + “ + •> + - + “ + “ + “ + - + - - + - + - + “ + - + - + - + 
C+-+” the sura of raagnitudes is |RD|**2 Difference is iRD-|**2 +-+ 
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DO 200 K=1-512 
DO 201 1=1,N 

OUTRCK,I)=ABS(0UTR(K,I)) 
OUTI(K,I)=ABS(OOTI(K,I)) 
0UTR(K,I)=0UTR(K,I) + 0UTI(K,I) 

C 

IF (0UTR(K,I) .GT, AMAX) THEN 
AMAX=OUTR(K,I) 

ENDIF 

IF (0UTR(K,I) .LT, AMIN) THEN 
AMIN=OUTR(K,I) 

ENDIF 

201 CONTINUE 
200 CONTINUE 
C+ 


PLOTTING . 

IF (C0NTR.EQ. ’C’)THEN 
GO TO 50 
ENDIF 

C CALL HEIGHTC0.28) 

CALL BSHIFT ( -0.2,-.25) 

CALL AREA2D(8,9) 

CALL VOLM3D{10,10,8) 

CALL HEADIN(TTL,100,1. ,3) 

CALL HEADINCSIGNAL,100,1.,3) 

CALL HEADINC TSMTH,100,1.,3) 

CALL ME3SAG(WTYPE,100,2.5,9. 3) 

CALL INTNOCZ ,'ABUT’,‘ABUT') 

CALL MESSAGC ' POINT'S)$’, 100, ‘ ABUT', ’ ABUT') 

CALL X3NAME('FREQUENCY AXIS$',100) 

CALL Y3NAME(’TIME AXIS$',100) 

CALL Z3NAME(' $',100) 

CALL VUANGL(-65,70,700) 

CALL XNONUM 
C CALL ZNONUM 

CALL MXIALFC ’ STANDARD','//') 

CALL MX2ALF('L/CGREEK','+') 

CALL ANGLE(-25.0) 

IF ( MODE .EQ. 2 ) THEN 
CALL MESSAGC +-P// ’,6,0. ,2. 3) 

ELSE 

CALL MESSAGC +0// ',5,0. ,2.3) 

ENDIF 

CALL ANGLE(-25.0) 

CALL MESSAGC +P# ’,5,4.9,0.15) 

CALL GRAF3D(-256,256,256,1,N,N,1.0*AMIN,0.5*(AMAX-AMIN), 
+ 1.0*AMAX) 

CALL SURMAT(OUTR,512,512,1,N,0. ) 

CALL ENDPL(O) 

CONTINUE 

50 IF r CONTR. NE. 'A')THEN 

CALL CONTORCOUTR,NX,N,TTL,SIGNAL,WTYPE,TSMTH,Z) 

ENDIF 
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WRITE(61,603)AMAX,AMIN 
603 FORMAT (IX,'MAXIMUM AMPLITUDE=',E14. 7 
+ /IX,'MINIMUM AMPLITUDE=',E14. 7) 

111 CONTINUE 

CALL DONEPL 
99 STOP 
END 

ssssssssssssssssssssssssssssssssss 

SUBROUTINE CONTOR(A,NX,NY,TITLE,SIGNL,V»NDW.TAVG,WLEN) 

THIS SUBROUTINE CONTOURS AN NX BY NY ARRAY OF REGULARLY SPACED POINTS. 
NOTE: THE ARRAY MUST BE REAL*4. 

A : SINGLE PRECISION NX BY NY ARRAY OF REGULARLY SPACED POINTS 
NX: NUMBER OF POINTS IN THE X-DIRECTION 
NY: NUMBER OF POINTS IN THE Y-DIRECTION 
ZINC: CONTOUR INTERVAL 

DIMENSION A(NX,NY) 

COMMON W0RK(50000) 

SET PARAMETERS FOR AXES: 

XORIG=-256. 

XSTP=256. 

XMAX=256. 

YORIG=0. 

YSTP=NY 

YMAX=NY 

SET CONTOUR LEVEL 
ZINC=AMAX/10 

CALL SETCLR('CYAN') 

SET PAGE AND PLOT SIZES, SET UP AXES FOR PLOT: 

CALL PAGE(8.5,11. 0) 

CALL BCOMON(50000) 

CALL AREA2D(6. 0,8. 0) 

LABEL AXES: 

CALL XNAMEC'FREQUENCY - AXIS $',100) 

CALL YNAMEC’TIME - AXIS $',100) 

CALL XNONUM 

CALL GRAF(XORIG,XSTP,XMAX,YORIG,YSTP,YMAX) 

CALL FRAME 

TITLE: 

CALL HEADINCCONTOUR PL0T$',100,1. ,4) 

CALL HEADINCTITLE,100,1.,4) 

CALL HEADINCSIGNL,100,1. ,4) 

CALL HEADINCTAVG,100,1.,4) 

CALL ANGLEC0.0) • ' 

CALL MESSAG(WNDW,100,1.5,-.7) 

CALL INTNOCWLEN ,'ABUT*,'ABUT') 
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CALL MESSAGC'POINTS)$*,100,’ABUT' / ABUT*) 

C 

C MAKE CONTOURS AND DRAW: 

CALL SETCLR('RED') 

CALL CONMINO.O) 

CALL CONANG(60. ) 

CALL CONLINC 0,'MYCON',’NOLABELS’,2,10) 

CALL CONMAK(A,NX,NY,’SCALE') 

CALL CONTURC1,'LABELS',’DRAW’) 

END PLOT: 

CALL ENDPL(O) 

RETURN 
END 

ssssssssssssssssssssssssssssssssss 

SUBROUTINE MYCONCRARAY,lARAY) 

DIMENSION RARAY(2),IARAY(1) 

THIS ROUTINE MAKES NEGATIVE CONTOURS DASHED AND THE ZERO LINE HEAVIER. 

CALL RESET('DASH') 

IF (RARAY(l) .GE. 0.) GO TO 10 
CALL DASH 
10 RARAY(2) = 1. 
lARAY(l) = 1 

IF (RARAY(l) .EQ. 0.) lARAY(l) = 2 
RETURN 
END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 


* * 

* CALL FFT(N,XTEMP,X) * 

* * 

* X - OUTPUT COMPLEX ARRAY CONTAINING FFT (1024) * 

* N - NUMBER OF POINTS * 

* XTMP - COMPLEX ARRAY CONTAINING DATA SAMPLES * 

* (starting at l,up to 1024) * 


*Vf*******yc***Vfyf**VfA***y«***vc*Vf*****Vf**A**yf**yf****rfr;v**** 


SUBROUTINE FFT(N,XTMP,X) 

COMPLEX X(1024),XTMP(1024),WTFAC.TMP 
M=INT(LOG10(FLOAT(N))/LOG10(2.)+0.5) 

EN = N 

PI = 4.0*ATAN(1. 0) 

DO 10 K=0,N-1 
NEWADR = 0 
MADDR = K 
DO 20 1=0,M-1 

LRMNDR = M0D(MADDR,2) 

NEWADR = NEWADR + LRMNDR*2**(M-1-I) 
MADDR = MADDR/2 
20 CONTINUE 


FFT00130 

FFT00140 

FFT00210 

FFT00270 

FFT00320 

FFT00330 

FFT00340 

FFT00350 

FFT00360 

FFT00370 

FFT00380 

FFT00390 
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X(NEWADR+1) = XTMP(K+1) 

10 CONTINUE 
DO 50 L=1,M 

ISPACE = 2**L 
S = N/ISPACE 
IWIDTH = ISPACE/2 
DO 40 J=0,(IWIDTH-1) 

R =: S*J 

ALPHA = 2,*PI*R/EN 

WTFAC = CMPLXC COS(ALPHA), -SIN(ALPHA)) 
DO 30 IT0P=J,N-2,ISPACE 
IBOT = ITOP + IWIDTH 
TMP = X(IB0T+1)*WTFAC 
X(IB0T+1) = X(IT0P+1) - TMP 
XCTOP+l) = X(IT0P+1) + TMP 
30 CONTINUE 

40 CONTINUE 
50 CONTINUE 
RETURN 
END 


FFT00400 

FFT00410 
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APPENDIX B. CROSS IPS 


Up to this point only autospectra have been discussed. Analysis of cross spectral 
characteristics of nonstationary phenomena can provide valuable information about the 
process. Equivalent to (1), a cross power spectral density is defined 


oo 


(63) 


For the case when x, and Xj are uncorrelated then 


(64) 


where n, is the mean value. If the data is correlated the energy resulting from cross 
spectral analysis can be complex. By examining Parseval's theorem in a more general 
context, 


r xi!)x]{t)dt=r 

= «V;(0), 

where /?x,,^(0) is not necessarily real nor is the cross correlation function (CCF) neces¬ 
sarily conjugate symmetric about /?,,,^(0). [Ref. 1] 

All the spectral estimators previously discussed are applicable if the ACF estimate 
is replaced by with the CCF estimate. The bias for cross spectra may be much larger 
than an equivalent autospectra where the point of maximum overlap occurs at lag zero. 
In practice, the location of maximum overlap is unknown for CCF's. One interpretation 
of cross PSD are of importance, the case where x, and Xj are two channels of a multi¬ 
channel system. The iP,^,^(/)l contains information concerning relative amplitudes at 
specific frequencies where the contains information concerning the lead or lag in 
phase between the tw'o channels. (Ref. 3] 
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Applying IPS in a multi-sensor environment leads to the following defining 
equation, 








( 66 ) 


where signal y could be a delayed, noisy version of signal x. Using a rapidly changing, 
linearly chirped pulse cM'mand implementing a windowed version of /PS,^, 
two sample cross spectra are shown in Figure 50 and Figure 51. 

The first spectra was created by beating the pulsed chirp x(r) against a delayed ver¬ 
sion of itself. The cross spectrum is shown in Figure 50 (a). In this case the delay is 
18 samples for a pulse 63 samples in duration. The cross spectrum can be seen in 
Figure 50 (b) and (c), where the spectral ridge over the interval corresponding to the 
absolute time of overlap. The maximum amplitude achieved on the cross-spectral sur¬ 
face is nearly the same as for the autospectra. The minimum however, is approximately 
44% greater in magnitude than that found on the corresponding auto-spectral surface. 
/PS.j, does not appear to provide information which can estimate delay in reception for 
this class of signal. 

The second cross spectrum considered examines the ability of IPS,^ to indicate cor¬ 
relation between a pulsed chirp and a Doppler-shifted version of itself. In this case, 

A'(r) = 'm lu ^0 

/IS / 2\ 


representing a shift of 50%. Figure 51 (c) shows an overlay of the two pulses. 
Figure 51 (b) is the contour for the cross spectrum. The peaks of the characteristically 
modulated ridge correspond to the line of overlap shown in (c). It is not clear if the cross 
spectrum of a linear chirp with a Doppler shifted version of itself yields information 
concerning the degree of coherence. The cross spectrum in Figure 51 could easily be 
interpreted as an autospectrum in which two, closelv-spaced parallel chirps are present. 
This initial investigation into the behavior of IPS,^ suggests that a more detailed exam¬ 
ination of its behavior is in order. 
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(c) IPS,^ amplitude plot 

Figure 50. Cross spectral analysis of a pulsed linear chirp 
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(c) IPS, overlayed contour plots of the original and Doppler-shifted pulse 


Figure 51. Cross spectral analysis for a Doppler-shifted linear pulse 







APPENDIX C. CUMULANT 


Assuming the received signal is corrupted by zero mean, Gaussian noise, examining 
the third-order moment or cumulant may yield information about the signal while sup¬ 
pressing the contributions of the noise. This potential processing gain is realized because 
the odd moments of a Gaussian process are identically zero. One way to implement IPS 
as an estimator of the cumulant is 


r»oo 

C«m./?S,Ro,(/;0 = y (U(/)|2.r(r-T) -b |/(/)|'4f + T)) 




Looking only at the effect on noiseless signal data, some preliminary results can be 
seen in Figure 52. All signals are analytic, therefor the magnitude squared term is al¬ 
ways unity. A comparison of the treatment of /PS,, can be made by referring to Chapter 
IV, Section C (3); Test Case Results. This method of processing the cumulant of an 
analytic signal does not appear to provide any useful results. Further research should 
examine the behavior using real signals or possibly forming the approximation of the 
cumulant in a different fashion. 



Figure 52. The cumulant of various analytic signals 
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